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1.  Introduction 

This  report  reviews  the  work  of  the  principal  investigator  as  published  in 
the  following  scientific  papers: 

1.  1.  Fried,  Nonlinear  finite  element  analysis  of  the  thin  shell  of 
revolution.  Submitted  for  publication  to  CMAME  (1983). 

2.  I.  Fried,  Orthogonal  trajectory  accession  to  the  nonlinear  equil¬ 
ibrium  curve.  Submitted  for  publication  to  CMAME  (1983). 

3.  1.  Fried,  On  unconditionally  stable  implicit  time  Integration 
methods  in  elastodynamlcs  and  heat  transfer.  Submitted  for  pub¬ 
lication  to  CMAME  (1983). 

4.  I.  Fried,  Reflections  on  the  computational  approximation  of  elastic 
incompressibility.  Computers  &  Structures  17,  161-168  (1983). 

5.  I.  Fried,  Nonlinear  finite  element  computation  of  the  equilibrium 
stability  and  motion  of  the  extensional  beam  and  ring.  CMAME  38, 

29-44  (1983). 

6.  I.  Fried,  Finite  element  computation  of  large  elastic  deformations. 
Proceedings  of  the  Bz  uiel  Conference  on  the  Mathematics  of  Finite 
Elements  and  Applications,  MEFLAP  IV,  J.R.  Whiteman,  Editor,  Academic 
Press,  143-159  (1982). 

7.  I.  Fried,  Finite  element  computation  of  large  rubber  membrane  deform¬ 
ations.  IJNME  18,  653-660  (1982). 

8.  1.  Fried,  Large  deformation  static  and  dynamic  finite  element  analysis 
of  extensible  cables.  Computers  &  Structures  15,  315-319  (1982). 

9.  I.  Fried,  Stability  and  equilibrium  of  the  straight  and  curved  elas- 
tica-finite  element  computation.  CMAME  28,  49-61  (1981). 


-vlt 


-2- 


10.  1.  Fried,  Nonlinear  finite  element  computation  of  the  equilibrium 
and  stability  of  the  circular  plate.  IJNME  17,  1427-1440  (1981). 

11.  I.  Fried,  Meaningful  existance  of  finite  element  solutions  to  off- 
limit  problems.  CMAME  22,  229-240  (1980). 

12.  1.  Fried,  Irregular  finite  element  meshes  in  elastodynamics . 

IJNME  15,  626-628  (1980). 

13.  I.  Fried,  On  the  optlonality  of  the  pointwise  accuracy  of  the  fin¬ 
ite  element  solution.  IJNME  15,  451-476  (1980). 

14.  I.  Fried,  Accuracy  of  string  element  mass  matrix.  CMAME  20, 

317-321  (1979). 

15.  I.  Fried  and  J.  Metzler,  SOR  vs.  conjugate  gradients  in  a  finite 
element  discretization.  IJNME  12,  1329-1342  (1978). 

16.  I.  Fried  and  J.  Metzler,  Conjugate  gradient  solution  of  a  finite 
element  elastic  problem  with  high  Poisson  ratio.  CMAME  15,  83-84 
(1978). 

17.  I.  Fried  and  J.  Metzler,  Displacement,  strain  and  stress  error 
nodal  lines  in  finite  elements.  Computers  &  Structures  9,  335-339 
(1978). 

18.  I.  Fried  and  D.S.  Malkus,  Energy  error  in  the  elastic  solution  when 
an  Incompressible  solid  is  assumed  compressible.  In  Formulation 
and  Computational  Algorithms,  K.J.  Bathe  et  al . ,  Editors,  MIT  Press 


131-139  (1976). 
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2.  Overview 

This  section  reviews  the  papers  listed  in  Section  1  in  ascending  order: 

1.  Energy  error...:  The  computational  problems  arising  in  the  finite  element 
simulation  of  incompressibility  have  occupied  the  author's  attention  since  1975. 
The  purpose  of  this  paper  is  to  show  that,  at  least  energetically,  the  exact 
solution  to  the  linear  three  dimensional  elastic  problem  depends  continuously 

on  Poisson's  ratio  v  .  This  means  that  taking  v  close  to  one  half  theoretical¬ 
ly  guarantees  a  close  analytic  solutions  to  the  incompressible  state.  This  is 
also  the  basis  for  the  residual  energy  balancing  technique.  Computational  dif¬ 
ficulties  concerning  convergence  and  conditioning  are  discussed  in  entry  4.  of 
Sec.  1. 

2.  Displacement,  strain  and  stress  error  nodal  lines...:  The  paper  shows 
computationally  the  existence  of  nodal  lines  in  two  dimensional  finite  elements 
on  which  the  error  between  the  computed  solution  and  the  exact  one  is  zero. 

3.  Conjugate  gradient  solution. . . :  Iterative  methods  for  the  solution  of 
the  large  stiffness  equation  set  up  with  finite  elements  has  some  distinct  ad¬ 
vantages  (and  admittedly  disadvantages)  over  direct  methods.  Iterative  methods 
do  not  require  the  explicit  assembly  of  the  global  stiffness  matrix  and  can 
operate  on  the  element  data  and  connectivity  Information  only  and  are  therefore 
immune  to  the  numbering  of  the  nodes  and  the  ordering  of  the  elements.  They  may 
require  the  minimum  of  data,  as  only  one  element  if  the  elements  repeat,  and 
naturally  take  advantage  of  symmetries  and  repeating  eigenvalues. 

The  paper  describes  the  application  of  the  conjugate  gradients  method  to 
the  solution  of  three  dimensional  elastic  problem  with  Poisson's  number  that  may 
be  close  to  0.5.  Scaling  of  the  global  stiffness  matrix  is  shown  to  be  highly 
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4.  SOR  vs.  conjugate  gradients. . . :  The  method  of  successive  over  relax¬ 
ation  has  achieved  great  popularity  in  mathematical  circles  for  the  Iterative 
solution  of  linear  systems  of  equations.  Its  main  drawback  is  that  it  requires 
an  elusive  factor  for  its  success.  It  is  shown  on  a  two  dimensional  heat  con¬ 
duction  problem  discretized  with  finite  elements  that  even  with  the  optimal  fac¬ 
tor  u  SOR  does  not  do  appreciably  better  than  the  conjugate  gradient  method. 

A  slight  change  in  the  optimal  «  causes  SOR  to  lose  to  CG  by  a  wide  margin. 

5.  Accuracy  of  string  element. . . :  Various  mass  distributions  are  dis¬ 
cussed  for  the  string  element  mass  matrix. 

6.  On  the  optimality...:  The  finite  element  method  produces  solutions 
that  are  energy  optimal.  Hie  paper  shows  that  the  theoretical  predictions  for 
the  pointwlse  optimality  of  the  finite  element  solution  are  correct.  Essen¬ 
tially  in  two  dimensional  second  order  problems  the  asymptotic  displacement 
error  for  linear  elements  is 

max  I u  -  ul  -  ch^in  ^ 

where  u  is  the  exact  solution  u  the  computed,  h  the  element  diameter,  and 
c  a  constant. 

The  paper  argues  that  the  asymptotic  error  estimate  is  impractical  since 
when  h  is  large  c  is  still  a  function  of  it  and  c(h)Jtn  does  not  behave 
like  in  ^  .  One  must  go  up  to  thousands  of  elements  to  be  able  to  numerically 
discover  the  in(^)  in  the  error  estimate. 

7.  Irregular  finite  element  mashes...:  The  paper  discusses  iom  basic 
problems  in  optimal  arrangement  of  finite  element  meshes  in  the  simple  context 
of  the  string  problem.  Setting  up  a  truly  optimal  mesh  of  finite  elements  is  a 
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very  difficult  task  even  in  the  simplest  of  problems  and  all  one  can  hope  is 
for  an  adaptive,  or  iterative,  procedure  for  the  mesh  Improvement.  The  paper 
discusses  the  mass  matrix  only  and  makes  the  following  interesting  conclusions: 

a.  The  optimal  mass  matrix  is  optimal  only  for  a  uniform  mesh.  Small 
deviation  from  uniformity  causes  drastic  accuracy  losses  with  this  element. 

b.  With  first  order  elements  the  grading  of  the  mesh  is  in  opposite  di¬ 
rections  for  the  lumped  and  the  consistent  mass  matrices.  The  consistent  mass 
matrix  is  also  insensitive  to  mesh  grading. 

c.  With  quadratic  elements  and  a  consistent  mass  matrix  the  mesh  grading 
is  opposite  to  that  for  the  linear  element,  and  the  problem  becomes  more  sen¬ 
sitive  to  mesh  variation. 

8.  Meaningful  existance. . . :  Considered  are  problems  that  are  theoreti¬ 
cally  off  limits  for  the  finite  element  method  such  as  problems  with  an  unlimi¬ 
ted  energy,  with  a  discontinuous  solution,  and  with  redundant  boundary  conditions. 
It  is  argued  that  the  finite  element  solution  in  all  these  cases  is  still  useful. 

9.  Nonlinear. . . circular  plate. . . :  Gauss  Integration  of  the  nonquadratic 
total  potential  energy  of  the  circular  plate  is  used  to  derive  the  element  tan¬ 
gent  stiffness  matrix  and  element  gradient  for  the  large  deformation  analysis  of 
the  plate.  Computations  are  made  for  the  plate  under  lateral  load  and  post  crit¬ 
ical  rim  thrust. 

10.  Stability  and  equilibrium. .. curved  elastics. . . :  A  slope  formulation 
of  the  elastics  is  dlscretised  with  quadratic  finite  elements.  The  element  tan¬ 
gent  stiffness  matrix  and  element  gradient  are  written  in  a  form  that  is  easily 
programmable  for  use  with  the  Newton-Bsphson  solution  of  the  nonlinear  stiffness 
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equatlon.  Load  and  stiffness  correction  solution  methods  are  also  considered  and 
the  importance  of  averaging  the  Iterates  is  demonstrated.  Very  large  deforma¬ 
tions  of  straight  and  bent  elastlcas  are  computed. 

11.  Large  deformation. . cables . . . :  To  extend  the  possible  boundary  conditions 
and  to  be  able  to  include  inertia  loads  the  cable  is  assumed  extensible.  A  quad¬ 
ratic-quadratic  large  displacement  element  is  formulated  through  discreet  sampling 
of  the  total  potential  energy.  Various  static  and  dynamic  cable  problems  are 
solved  and  the  difficulties  in  computing  the  tension  from  the  constitutive  equa¬ 
tion  is  pointed  out.  An  experiment  with  a  falling  cable  shows  excellent  agree¬ 
ment  with  the  computer  model. 

12.  Finite  element. . .rubber  membrane. . . :  Numerical  sailing  of  the  total 
potential  energy  is  used  to  derive  the  quadratic  element  tangent  stiffness  and 
element  gradient  for  the  largely  deformed  axisysnetrlcal ,  Mooney,  rubber  membrane. 
Computations  are  made  for  the  stretched  and  Inflated  disc,  the  inflated  torus, 
and  the  bulging  tube  under  Internal  pressure.  The  convergence  of  the  Newton- 
Raphson  method  near  a  critical  point  is  discussed. 

13.  Finite  element ... large  elastic  deformation. . . :  Creation  of  the  element 
tangent  stiffness  matrix  and  element  gradient  for  the  computation  of  large  elastic 
deformation  is  discussed.  The  procedure  does  not  require  the  elements  to  be  small 
in  the  sense  of  approximating  an  arc,  and  is  otherwise  analogous  to  the  linear 
finite  element  method  except  that  the  Newton-Raphson  (or  its  modification)  method 
must  be  called  to  solve  the  nonlinear  stiffness  equation. 

14.  Nonlinear. .equilibrium. . . extenaional  beam. . . :  The  element  data  for  the 
large  deformation  analysis  of  the  curved  extensible  rod  are  given.  Numerical 
tests  are  done  with  the  element  to  compute  the  large  deflection  of  a  cantilever. 
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the  pressing  of  a  ring  by  polar  forces,  the  deformation  of  a  circular  ring  by  ex¬ 
ternal  pressure;  the  motion  of  simply  supported  and  free  beams,  the  vibration  of 
the  ring,  and  the  bending  of  a  cantilever  by  a  follower  force. 

15.  Reflections  on... elastic  incompressibility...:  Introduction  of  in¬ 
compressibility  into  the  finite  element  model  is  still  of  great  current  interest. 
It  is  argued  that  the  gradual  Increase  of  the  bulk  modulus  coupled  with  a  mesh 
refinement  is  the  most  sensible  thing  to  do  both  from  the  theoretical  as  well  as 
the  practical  point  of  view. 

16.  On  Unconditionally  stable. . . :  Unconditionally  stable  explicit  time  in¬ 

tegration  methods  are  of  the  greatest  interest  in  the  finite  element  analysis  of 
elastodynamics  and  heat  transfer.  Several  such  algorithms  have  recently  been 
published.  It  is  shown  that  the  time  and  space  errors  become  coupled  and  that 
the  space  mesh  reduction  may  cause  a  disasterous  loss  of  accuracy. 

17.  Orthogonal  trajectory. .. :  The  Rlks-Wempner  method  for  correcting  the 
equilibrium  equation  is  modified  to  remove  the  need  for  an  explicit  load-dis¬ 
placement  constraint. 

18.  Nonlinear. . .shell  of  revolution. . . :  A  cubic-cubic  element  is  developed 
for  the  large  deformation  analysis  of  the  axisymmetrical  shell  of  revolution. 
Explicit  formulas  are  given  for  the  element  tangent  stiffness  matrix  and  element 
gradient.  The  bending  of  a  thin  spherical  shell  under  the  combined  action  of 
polar  forces  and  a  surface  pressure  is  computed.  The  high  computational  price 
of  the  Newton-Raphson  method  is  noted. 
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Integration  by  parte  and  addition  of  the  resulting  form  of  equation  (19) 
bade  to  %(U)  results  in  the  cancellation  of  the  last  two  Integrals  In 
equation  (18)  and  ♦(U)  becomes  now 
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DISPLACEMENTS,  STRAIN  AND  STRESS  ERROR  NODAL 
LINES  IN  FINITE  ELEMENTSt 

Isaac  Pried 

Department  of  Mtikmiin.  Romo*  University.  Boston.  MA  02215.  U.S.A. 
aad 

J.  A.  Metzler 

.  Department  of  Malhcmalics.  Draw  University,  Madison.  NJ  07940.  USA. 

(XfcnW  II  September  1977;  reeehet  for  pabtkuioo  I  December  1977) 

Abstract— Numerical  experiments  show  dial  Ike  error  ia  ike  computed  fcmte  element  solution  and  in  faimiw 
vanish  on  typical  error  nodal  lines  inside  each  Inite  element.  A  ibeoreiical  explanation  it  five*  to  this  phenomenon. 
wfcich  was  previously  discovered  for  distinct,  special  points.  Systematic  dassiScatioo  of  tkese  lines  far  diErirnt 
dement  types  and  problems  appears  to  be  a  wonky  undertaking. 
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>  Since  the  oricinal  observation  by  Bariow(l)  that  the 
.  (mate  element  stresses  computed  in  rectangular  elements 
i  d  the  Gauss  points  are  superior  in  accuracy,  there  have 
i  keen  repons  on  different  such  special  points  in 
iectangks(2)  and  triangles(3].  Using  the  best  energy  flt 
.  technique  [4],  whereby  minimization  of  the  error  in  the 
;  energy  b  carried  out  over  a  single  typical  clement  with 
I  in  assumed  polynomial  exact  solution  that  assures  a 
(  tfobally  admissible  finite  element  solution,  it  is  shown 
here  that  there  arc  entire  lines  (surfaces  in  space)  inside 
dtt  element  on  which  accuracy  b  superior  and  hence  the 
!  nriety  of  observed  special  points. 

These  nodal  lines  depend  on  the  dement,  original 
4  fmbfcm  and  eventually  location,  as  shown  theoretically 
;  Md  numerically  in  this  paper. 


•  momn-s  EQUATION 

Here  the  finite  element  solution  i  b  obtained  from  the 
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:  "equivalently 

|  w(l)- +(«,-!,)*]  dxdy  (2) 

^  u  and  4  are  the  exact  solution  and  the  finite 
.  ‘•■tat  trial  function,  respectivdy. 
f  for  tic  bthnear.  reactaagular  finite  element  4 
■  **+eiy+usxy.  -I  *x.y  *  I,  inside  each  ekment  and 
’'{"■tiufily  Ihe  choke  n  m  n*x* + *,xy + mij*  b  made, 
"■■imizaiion  of  (2)  yidds  a,  -  0, «,  ■  0  and  «,  -  o,  such 
7*  d»n»+n,ey.  Equating  a,  with  d„  I, 
"e*Ms  Ma0,  and  ■  #,  and  dm  arms  It,* 
^*lttndh,»i,-l,ui  changing  signs  on  the  lines 
*"•  and  y-0,  respectively.  Choosing  dm  hoe 
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parameter  a««0  makes  the  finite  dement  solution  the 
interpolate  to  the  quadratic,  meaning  that  it  b  extendable 
over  neighboring  elements.  Also,  since  any  smooth  solu¬ 
tion  is.  by  Taylor's  theorem,  almost  a  polynomial  inside 
each  ekment,  it  b  expected  of  the  errors  to  change  signs 
on  x  *0  and  y  «  0  in  general.  One  b  also  tempted 
conjecture  (as  was  done  m  [4])  that  the  best  place 
compute  S  b  at  the  nodes.  Numerical  experiments, 
however,  refute  thb.  at  bast  in  the  biinear  cate. 

For  the  biquadratic  dement  fi-na+ngr  +  nxy  +  na* 
♦  n«xy  +  e»y*+na*y  +  nay*  ♦  es*'y\  -I*x,y*  |. 
and  consequently  the  choice  rrwanx’+nix’y  +  ntzy* 
+  •»»’  is  made.  Similar  argumnets  bad  hen  to  the  con¬ 
clusion  that  tn.  changes  signs  along  the  lines  n  ■  ±V<3/3> 
and  tn,  along  y  -  ±V(3/3)  For  Sn  one  b  again  led  to 
expect  the  best  accuracy  at  the  nodes. 

To  numerically  observe  all  thb  and  to  toe  the  effects 
-bf-aumerical  integration. the  aquation  Um*e„»  I  nfib 
solved  in  a  square  region  with  a  on  its  boundary. 
Finite  element  calculations  ate  carried  out  with  biinear 
and  biquadratic  dements  (because  of  symmetry  only  one 
quarter  b  considered.)  The  error  nodal  lines  of  fin.  and 
In,  for  a  7x7  mesh  of  bilinear  dements  and  the  shaded 
ekment  as  in  Fig.  1(a),  are  shown  m  Figs.  1(b)  and  (ck 
Figures  1(d)  and  (e)  show  the  nodal  lines  of  fn  aad  la, 
for  the  biquadratic  ekment  (3x3  mesh  aad  an  ekment  at 
the  middb  of  the  quadrant),  with  exact  integration  and  a 
2x2  (broken  line)  Gauss  integration  scheme.  If  d.  or  I, 
are  confuted  at  the  Gauss  point  then  the  2x2  in¬ 
tegration  results,  as  can  be  teen  tram  Tabb  I,  ia  a 
certain  loss  of  accuracy. 

Here  one  b  caled  to  minimize  (Poisson's  ratio  »>•).' 
The  potential  energy 
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Introduction 

The  conjugate  gradients  method  has  been  shown  [11  to  be  an  attractive  technique  for  the  solu¬ 
tion  of  the  linear  algebraic  system  that  is  produced  by  finite  elements.  It  has  been  noticed  however 
[21  that  convergence  depends  on  the  eigenvalue  spread  in  the  global  stiffness  matrix,  and  that  an 
ill-conditioned  system  leads  to  poor  convergence.  Scaling  has  been  shown  to  improve  the  per¬ 
formance  of  the  conjugate  gradient  method  in  plate  problems.  In  the  present  paper  this  method 
is  applied  to  an  elastic  problem  to  study  its  functioning  in  the  presence  of  high  Poisson  ratios  that 
cause  a  deterioration  [3]  in  the  condition  of  the  global  stiffness  matrix.  It  is  also  shown  here  that 
scaling  has  a  considerable  beneficial  effect  on  the  convergence  of  conjugate  gradients  also  in  this 
case. 


I.  The  elastic  problem  and  its  discretization 

We  propose  to  use  conjugate  gradients  to  solve  the  problem  of  a  rotating  hollow  elastic  sphere 
[4]  discretized  with  finite  elements.  Because  of  symmetry  only  a  quarter  of  the  sphere  need  be 
considered,  and  we  divide  the  arc  and  radius  into  Ne  equal  parts  to  form  a  mesh  with  N?  “aquare” 
finite  elements.  A  biquadratic  interpolation  is  adopted  over  each  element,  which  is  thus  associated 
with  nine  nodal  points. 


2.  Numerical  computations 

Actual  numerical  computations  were  carried  out  with  Ne  =  4,  5, 6  and  7  elements  per  side,  once 

*  Work  supported  by  ONR  contract  No.  ONR-N00014-76-C-36. 
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without  scaling  and  once  with  the  global  stiffness  matrix  K  symmetrically  scaled  so  that  Ku  =  1 . 
The  algorithm  was  terminated  in  each  case  when  the  change  in  the  quadratic  form  tt(jc)  =  \  xxKx 
-x*/ which  is  minimized  reached  the  machine  accuracy.  The  results  of  these  numerical  computa¬ 
tions  are  listed  in  the  tables  below.  Notice  in  these  tables  the  substantial  savings  with  scaling.  The 
following  notation  will  be  used: 
v  Poisson’s  ratio 
Ne  number  of  elements  per  side 
N  size  of  the  linear  system 
Nit  number  of  iterations  required  for  convergence 
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IMRODI  runs 

hvcti  though  finite  elements  arc  being  accepted  now  .is  ills-  most  appiopnatc  w a\  to  discrcti/c 
problems  which  tan  bt  formulattil  v  ai  i  utionallv .  the  question  ol  how  to  .solve  the  resulting 
algebraic  svstents  is  \et  ninth  in  debate .  I  hi-  l.ucc  coinpulatiun.il  packages  have  almost 
exclusivclv  adopted  direct  nicthods  ol  solution  I'asnl  on  Gauss  diminution  and  triangular 
factoi i/ation  Ihe  reasons  lot  this  choice  ate  mam  the  need  to  solve  the  same  svstem  with 
(Itlleret’t  t  igltt  ha  ml  v  eclors.  the  direct  method  s  small  sensitiv  itv  1. 1  ton  ml  oil  errors  and  then 
termination  in  a  finite  niimhct  of  steps.  Direct  methods  sutler,  however,  from  the  disadvantage 
thill  they  operate  on  a  two- diittetision.il  ai  i  av  whu  h  in  the  ease  ol  finite  elements  ts  sparse  but 
with  no  distinct  spaiseness  pattern.  As  a  tesult  the  efhnt  to  make  direct  methods  etlieient 
becomes  a  considerable  progiummmu  exercise  to  avoid  zero  opoiations  that  depend  also  on 
the  numbering  ol  the  nodal  points 

Iterative  methods  are  subtlei  to  applv ;  one  has  to  know  when  to  stop  them  and  their 
ethcicncv  is  a  hidden  tunction  of  the  piopeilies  of  the  global  stillness  matrix.  On  the  other 
hand  their  programming  is  consuleiublv  simpler  than  that  ot  dit  .ut  methods  and  the  sparseness 
of  the  matrix  can  be  accounted  for  in  an  uncomplicated  vv.iv  since  thov  do  not  operate  on  a 
two-dimensional  arr.iv  but  lather  ivqunv  ;is  ;i  basic  opciation  onlv  the  multiplication  of  a 
vector  bv  it  niiitrix. 

Among  iterative  methods.  Successive  <  >v et iv'lax.ition  lias  icccivcd  extensive  coverage  in 
the  context  of  finite  difference  equations  lm  some  tcason  Conjugate  Gradients  never 
became  that  prominent  This  method,  hovvrvei  is  tei'huieallv  v.  iv  well  suited'  to  solve  Ihe 
tmile  element  ult'ibi.iii  equalious  Alllimii'li  biith  ol  these  m,  tli.nl.  Millet  I  tom  tiiinealioii  and 
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round-off  errors,  it  is  the  author's  experience  that  this  b  not  a  serious  prdMer 
system  is  particularly  ill-conditioned.  In  the  cases  where  tlr.s  :«  a  problem, 
recommend  scaling  of  the  type  described  in  References  3,  4.  »  ;  ^ 

In  this  paper  we  apply  both  SOR  and  CG  to  a  two-dimensional  stationary  heat  £wnlluctloo>> 
problem  disereli/ed  with  finite  elements  and  show  that  even  undei  optimal  conditions,  that  in 
fact  rarely  materialize.  SOR  converges  only  slightb  faster  than  the  more  conveniently  pro¬ 
grammed  CG. 


SUCCESSIVE  OVER  RELAXATION 

Let  us  recall  the  SOR  algorithm.  We  wish  to  solve  the  linear  algebraic  system 

Kx  f  (I) 

and  with  positive  definite  and  symmetric  K,  which  we  decompose  in  the  form 

K  I)  v  1.  *  I.'  (2) 

where  I)  is  (block)  diagonal  and  I.  the  remaining  strictly  lower  triangular  matrix.  We  write  now 

(I)  +  <ul.|X|  |  +(|  <u)l)|x,,  i  <ef  (3) 

in  which  to  is  the  overrelaxation  factor.  The  successive  overrelaxation  scheme  in  equation  (3) 
has  the  iteration  matrix 

f.,„  |1>  t  «»1.)  '|  wl.'  +(1  <•;)!>)  (4) 

and  the  basic  problem  of  SOR  is  to  locate  the  optimal  to  that  minimizes  the  spectral  radius 

(»(/.,.,)  of  the  matrix  No  sure  way  exists  to  find  the  optimal  »».  w  hich  is  the  greatest  flaw  of 
SOR  especially  in  view  of  the  fact  that  the  efficiency  of  the  method  strongly  depends  on  <o. 


CONJUGATE  GRADIENTS 
Here  the  algebraic  system  in  equation  (I )  is  solved  by 

r„-=f-Kx,,.  p„  r„ 
a,  =  p'r./p.'Kp, 
x, ,  i  =  x,  *-«,p, 
r,.i  r,  a,Kp, 

P,  r.M/r.'r, 

P..,  =r,.,  +  P.p, 

No  elusive  factor  is  needed  in  this  algorithm  and  the  only  matrix  vector  multiplication  that 
appears  in  it  is  Kp  How  to  carry  out  this  algorithm  with  finite  elements  is  discussed  in 
Reference  .3.  What  we  want  to  do  here  is  only  to  compare  the  two  algorithms  on  a  given 
realistic  problem. 


TEST  CASE 

Figure  I  shows  half  the  cross  section  of  a  double  barrelled  steam  pipe  of  which  abedtf  is  a  line 
of  symmetry.  The  Ittrger  bore  b<  carries  steam  at  a  constant  temperature  SfKfF  while  the 
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smaller  is  designed  Id  hose  a  ihermomeiei  .mil  we  assume  lire  aieli  he  insulated.  Heal  is  losi 
from  the  oulei  surface  til  of  the  manifold  hv  initiation.  I'he  etoss  section  is  divided  into  140 
triangular  elements  with  the  tcmpcinture  taken  to  he  linear  o\et  each  one  of  them.  This 
discretization  gives  iise  to  W  unknown  tcmpeiatuic  values  at  the  vertices  of  all  elements. 

I  hc  iteration  matiix  in  equation  (4  )  was  set  rip  and  its  sjvctial  tadms  ,i f  computed  as 
a  function  of  m  anil  plotted  in  I  iemc  2.  It  is  seen  that  the  optimal  <n.  i»  is  I  (W.  Next  the  heat 
conduction  piohlem  was  solver)  first  with  (  (  i  and  then  with  NOR  using  this  re,,,,!.  Trr  compare 
the  convergence  in  both  cases  the  iterated  temperature  at  point  I.  T.  is  plotted  in  Figure  4 
versus  the  number  ot  iterations  V  ,  When  m,,,„  is  useil  in  SOR  it  is  seen  to  he  slightlv  faster 
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Various  element  mass  matrices  are  considered  for  the  two-nodal-point  string  element.  From  finite 
difference  arguments  it  is  known  that  an  optimal  element  mass  matrix  exists  -  equal  to  half  the  sum  of 
the  lumped  and  consistent  mass  matrices  -  which  causes  an  error  of  only  0(h4)  in  the  computed  natural 
frequencies  instead  of  the  0(fi2)  that  is  obtained  with  the  consistent  matrix.  It  is  shown  here  that 
nonuniformity  in  the  mesh  destroys  this  optimality  and  has  also  an  adverse  effect  on  the  accuracy  of  the 
frequencies  computed  with  the  lumped  element  mass  matrix. 


1.  Uniform  mesh 

Let  us  write  the  element  stiffness  matrix  kt  and  element  mass  matrix  m,  for  the  two-nodal-point 
string  element  as 


where  the  element  size  h  may  vary  from  element  to  element.  When  a  =  2/3  and  a  =  1,  m,  in 
eq.  (1)  becomes  the  (variationally)  consistent  and  lumped  element  mass  matrix,  respectively. 
According  to  the  minimax  principle  [1]  all  the  natural  frequencies  computed  with  the 
consistent  m,  are  above  the  corresponding  exact  ones,  and  the  error  in  them  is  0 (h2). 

We  are  interested  in  studying  the  effect  of  a  on  the  accuracy  of  the  computed  natural 
frequencies  of  the  fixed  string  with  the  hope  of  discovering  an  optimal  a.  We  shall  do  it  by 
finding  a  closed  form  expression  for  the  approximate  natural  frequencies  of  the  fixed  string. 
Assembly  of  two  neighbouring  elements  produces  the  finite  difference  equation  for  the  /th 
node: 


-M/-,  +  2 My  -  u;>,  =  |  aV[(1  -  a)Uj~i  +  2 aUj  +  (1  -  a)u/+I],  (2) 

where  to  =  Va  is  the  approximate  natural  frequency,  where  j  =  1, 2, . . .  N,  and  h  =  1/(N  +  1). 
Eq.  (2)  is  solved  by  u,  =  cz',  which  when  substituted  in  eq.  (2)  leads  to  the  characteristic 
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equation  for  z: 

z2-2bz  +  1  =  0, 

L  1-an  1  2 

b  1  +  fi(l  —  a)'  *  2kh  ' 

(3) 

or 

2  =  6±ViM. 

(4) 

Since  the  eigenfunctions  of  the  fixed  string  are  trigonometric  functions,  we  expect  the  roots  of 
the  characteristic  eq.  (3)  to  be  complex.  The  condition  for  this  is  that  b2-  l  <0,  or 

0<m(2«-1)<2,  (5) 

assuming  that  2a  -  1  >  0.  We  know  [1]  that  the  eigenvalues  f  of  the  global  system  Ku  =  AAfu  are 
bounded  by  the  eigenvalues  of  the  element  system 


k'U,  =  A 


(6) 


such  that 


min  {A  t}  ^  A  <  max  {A  (7) 

e  e 

where  A'  and  A'  are  the  first  (lowest)  and  nth  (highest)  eigenvalues  of  eq.  (6).  Here,  with  A, 
and  me  given  in  eq.  (1)  A'  =  0  with  the  corresponding  element  eigenvector  u\  =  [1,  1J,  and 
A 5  —  4/[(2a  -  l)/i2]  with  the  corresponding  element  eigenvector  uL  =  [1,  -1].  Hence,  according 
to  eq.  (7),  O^/i  <21  (2a  -  1).  But  /i  =  2/(2a  -  1)  occurs  only  when  the  string  is  free-free  for 
the  global  eigenvector  u1  =  [1,  -1,  1,  -1, . . .].  Fixation  of  the  end  points  reduces  the  maximal 
natural  frequency  to  fi<  2/(2a  -  1),  and  z  in  eq.  (4)  is  indeed  complex.  Because  the  free  term 
in  the  characteristic  eq.  (3)  equals  1,  |z|  =  1,  and  we  may  write  the  complex  conjugate  roots 
of  this  equation  as  z  =  cos  0±isin  $.  Then  u,  =  Cie1'*  +  c2e  ,/*,  i  ->  or  w,  = 

Ci  cos  jd  +  c2  sin  jO,  j  =  0, 1, 2, . . .  N  +  1 .  The  end  condition  u0  =  0  is  satisfied  witti  Ci  =  0,  and  the 
end  condition  uN+ 1  =  0  with 


(N  +  1)0  =  7r,  2jt,  . . .  nir.  (8) 

The  yth  entry  in  the  nth  eigenvector  therefore  becomes 

M  =  sin^j-  (9) 

From  eq.  (4)  we  have  that  cos  8  =  b,  or 

_ _ 1  -  cos  0 

M  1  +  (a  -  1)0  -  cos  0)' 

i 


(10) 
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which  (with  the  identity  1  -  cos  8-2  sin3  \0)  simplifies  into 


u>h  = 


_ 2  sin  {d _ 

[1  +  2(a  -  l)sin3  \$]tn 


(ID 


Expansion  of  the  right-hand  side  of  eq.  (11)  in  terms  of  6  yields  the  following  error 
expression  for  a): 


^-l  =  -s(6a-5)03  +  O(04).  (12) 

ai 

When  a  =  4/6  (i.e.  when  me  is  consistent)  and  when  a  =  6/6  (i.e.  when  me  is  lumped),  the  error 
in  the  computed  frequencies  is  0 (h2).  But  when  a  =  5/6,  the  error  (according  to  eq.  (12)) 
decreases  to  0(/i4)  (see  [2],  [3]).  In  view  of  this  we  shall  term  m,  with  a  =  5/6  optimal.  It  is 
interesting  that  the  optimal  element  mass  matrix  is  obtained  from  the  physical  lumping 
method  [4]  by  considering  half  the  element  mass  uniformly  distributed  and  another  half 
lumped  at  the  ends. 

Fig.  1  compares  the  frequency  errors  for  the  different  choices  of  a  in  m,  in  eq.  (1)  for  the 
complete  spectrum  0  ^  £  =  \nnh  <  \n. 


2.  Nonuntform  mesh 

Accuracy  predictions  obtained  by  finite  difference  arguments  over  a  uniform  mesh  are 
susceptible  to  deterioration  when  the  mesh  is  not  uniform.  To  see  what  happens  to  the 
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accuracy  of  the  frequencies  of  the  fixed  string  computed  with  the  various  element  mass 
matrices  me,  we  numerically  solve  the  string  eigenproblem: 

u"  +  to2u  =  0,  0<x<l  (13) 

u(0)=  «'(j)  =  0 

for  the  fundamental  frequency  a>i. 

According  to  eq.  (13)  u "  is  proportional  to  u,  suggesting  a  nonuniform  mesh  with  elements 
near  x  =  0  larger  than  those  near  x  =  We  decide  to  grade  the  elements  according  to 


hj  —  e  sin' 


-1/2  . 


tt; 


2 (N  + 1)’ 


/=  1,2,  ...N  +  l, 


(14) 


where  e  is  fixed  by  the  condition  that  h,  +  h2  + - 1-  hN+]  =  j.  The  graded  mesh  thus  obtained 

is  shown  in  fig.  2  for  the  case  of  12  elements  (Ne  =  N  +  1  =  12). 


Fig.  2.  Convergence  of  the  fundamental  frequency  At  in  a  fixed  string  that  is  discretized  with  (a)  a  uniform  mesh 
( - )  and  (b)  a  nonuniform  mesh  ( - ). 


Computation  of  the  fundamental  frequency  to2  =  A,  is  now  carried  out  first  with  a  uniform 
mesh,  and  the  error  reduction  with  the  number  of  elements  Ne  is  indicated  in  fig.  2  by  a 
broken  line.  As  predicted  in  the  previous  section,  the  magnitudes  of  the  error  in  A,  computed 
with  the  lumped  and  consistent  mass  matrices  are  nearly  equal.  The  optimal  mt  produces  a  far 
more  accurate  Ai  with  a  faster  diminishing  error.  The  erratic  convergence  of  At  computed  with 
the  optimal  mt  (as  seen  in  fig.  2)  is  due  to  round-off  errors,  the  calculation  being  carried  out  on 
a  computer  with  about  8  significant  digits. 

The  nonuniform  mesh,  graded  according  to  eq.  (14),  produces  a  slightly  more  accurate  A1 
when  the  consistent  m,  is  used.  On  the  other  hand,  the  effect  of  nonuniformity  on  the  accuracy 
of  A,  computed  with  the  optimal  me  is  devastating.  Superconvergence  is  not  only  wiped  out, 
but  the  error  in  Ai  is  not  even  proportional  to  Ne~2  in  the  range  N,  12. 
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The  grading  rule  (14)  that  is  based  on  variational  error  estimates  produces  also  a  less  accurate  A*, 
with  the  lumped  mt.  It  is  observed,  however,  that  to  improve  Ai  with  the  lumped  m,  the  mesh 
should  be  with  smaller  elements  near  the  ends  than  in  the  middle. 
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ON  THE  OPTIMALITY  OF  THE  POINTWISE  ACCURACY  OF  THE 
FINITE  ELEMENT  SOLUTION! 

ISAAC  FRIED! 

Department  of  Aeronautical  Engineering.  Technion-Israet  Institute  of  Technology,  Haifa,  Israel 


SUMMARY 

Closed  form  finite  element  solutions  are  obtained  for  the  uniformly  loaded  membrane  in  R",  n  =  2,3 
discretized  with  first -order  elements.  It  is  verified  on  this  model  that  the  pointwise  error  in  the  displacement 
is  0( A2  log  1/h).  Slopes  converge  pointwise  at  the  optimal  rate  0(/i). 


INTRODUCTION 


Since  the  finite  element  solution  is  obtained  by  the  minimization  of  the  total  potential  energy, 
this  solution  is  optimal  in  the  energy  norm.  If  the  element  shape  functions  include  a  complete 
polynomial  of  order  p,  then  the  error  in  the  finite  element  solution  to  second-order  (membrane) 
problems  is  0(hp).  Energetically,  this  is  the  best  approximation  achievable  with  piecewise 
polynomials  of  order  p.  But  what  about  the  pointwise  error  in  the  finite  element  solution?  The 
best  pointwise  displacement  error  possible  with  elements  of  order  p  is  0(/i',+  1).  Recently,  Scott1 
showed  that  for  the  two-dimensional  membrane  the  pointwise  error  in  the  computed  solution  u 
is  actually 


0(h2  log  1/h), 

0(/ip+1), 


p=\ 

p&2 


(1) 


Nitsche2  showed  that  with  linear  element  the  pointwise  error  is  0(/t2  log  l/h)  also  in  three 
dimensions. 

Here  the  uniformly  loaded  membrane  problem  in  R",  n  -  2,  3 


-  r'~’'(rn~'u')'  =  2n,  0<r<l 

«'(0)  =  «(1)  =  0 


(2) 


that  is  solved  by  u  =  1  —  r2  is  considered  (see  also  Reference  3).  With  the  total  potential  energy 

*•(«)=[  (i«’*-2n«)r"-'dr  (3) 

Jo 


a  finite  element  solution  u,  consisting  of  piecewise  linear  elements,  is  obtained  from  equation  (1 ) 
in  a  closed  form.  It  is  verified  from  this  solution  that,  in  fact,  asymptotically  ||u  —  m||od  = 
0(h 2  logl /h).  Slopes  converge,  according  to  the  present  computational  model,  at  the  optimal 
rate  ||«’-  fi'll*  =  0(h),  both  in  two  and  three  dimensions. 
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CIRCULAR  MEMBRANE 


Discretization  is  performed  with  N  =  \/h  linear  elements.  With  n  =  2,  the  total  potential  energy 
in  equation  (3)  provides  the  element  stiffness  matrix  k,  and  element  load  vector  f,  as 


from  which  the  global  data 


(4) 


K  =  i 


'  1  -1 

'  f 

-1  4 

-3 

6 

-3 

8 

-5 

-5 

12  -7 

ll 

12 

18 

• 

(5) 


is  assembled.  Because  of  the  regular  structure  of  K  and  f  the  global  system  Ku  =  f  is  readily 
solved  to  yield 


A)3 


.r-u-ir 


»  =  1,2 . N 


but 


/•;-</- 1>1  3/2-3 /+!  ,/  1  \ 

/5-(/-l)2  2/-1  4l6/  3  2/  —  1  / 


(6) 


(7) 


and  consequently  equation  (6)  becomes 

(8) 

/-■  2/  - 1 

where  r,  »  h(i  - 1).  Figure  1  traces  the  error  distribution  u-u,  where  u  =  1  -  r1  and  u  is  linearly 
interpolated  inside  the  element  from  the  computed  nodal  values  u,  in  equation  (8),  over  a 
seven-element  mesh.  The  maximum  pointwise  error  occurs  at  r  =  0  (i.e.  i«l,  r,  =  0)  and 


Figure  i .  Error  distribution  u-A  over  ■  seven  linear  elements  discretisation  of  the  uniformly  loaded  circular  membrane 
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and  when  j  =  (N  +  i  )/2,  for  instance,  u'-u‘=  \N  2;  superconvergence  takes  place.  The  maxi¬ 
mum  error  in  u'  is  at  r  =  0,  and 

max  \u'-u'\  =  *h  (15) 

Slopes  converge  pointwise  at  the  optimal  rate  0 (h ). 


SPHERICAL  MEMBRANE 


With  n  —  3,  the  total  potential  energy  in  equation  (i)  yields  the  element  data  for  the  first-order 
element  as 


maths  grids 


k,  =  jh(3e2-3e  +  1)J_|  }]. 


6e2-8e  +  3 

6e2  — 4e  + 1 


] 


from  which  the  global  data 


1-1 

r 

-1  8  -7 

14 

-7  26  -19 

50 

-19  56  -37 

110 

.  f-W 

-37  98  -61 

124 

. 

. 

(16) 


(17) 


it  assembled. 

A  closed-form  solution  to  K&  ■*  /  is  readily  obtained  and  at  the  rth  node 


*,-i  h2 


N 


l 

i-i 


/*-(/- D4 
/*-</- 1># 


(18) 
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But 


/*-(/  - 1)4  4/3-6j2+4j-  1 
P-O'-i)3  3/ 2  -  3/  + 1 


2/  - 1 
3/2-3/+l 


) 


(19) 


and  consequently,  equation  (18)  becomes 


M.  = 


\-r]+W 


f  2/-1 

,t,3/2-3/  +  l 


(20) 


Figure  4  traces  the  error  w  -  u,  where  u  =  l-r1  and  where  u  is  linearly  interpolated  from  the 
nodal  values  u*  given  in  equation  (18),  over  a  seven-element  discretization.  The  largest 


Figure  4.  Sphere,  displacement  error  distribution  u  -  u 


pointwise  error  in  «  occurs  at  r  =  0  (i.e.  i  =  1),  and  according  to  equation  (20) 


llw  —  m||oo  =  2*2  I 

/-I  Jj 

(21) 

Again,  equation  (21)  is  written  as 

||«  -  m||«  =  ic(Af)Af’2  log  N 

(22) 

where  c(N)  is  graphed  in  Figure  2.  As  N->  oo,  c(N)-*  2/3,  but  so  slowly  as  to  make  the 
asymptotical  error  estimate  ||u  -  i3|U>  =  0 (N~2  log  N)  of  theoretical  interest  only. 

The  computed  slope  inside  each  element  is  here 

(23) 

and  the  error  distribution  u'  -  6'  is  graphed  in  Figure  5.  Change  of  sign  ofu'-u'  takes  place,  also 
here,  inside  each  element  and  the  error  nodal  points  for  u'-u'  are  at 
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The  maximum  of  \u'  -  u'\  occurs  at  r  =  0,  and 

max|w'-«'|  =  2A  (25) 

Also  in  three  dimensions  the  computed  slopes  converge  pointwise  at  the  optimal  rate  0(/i ). 
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IRRFGULAR  FINITE  ELEMENT  MESHES 
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SUMMARY 

Grading  rules  for  consistent  arid  lumped  elements  ire  opposite  Optimal  mass  matrices  on  finite  difference 
considerations  are  very  sensitive  to  mesh  grading.  Scnsitivitv  of  the  solution  accuracy  to  the  mesh  ratio 
increases  vvuh  element  urdet. 


MRS  I  ORDI  R  I  I  .FMI  S'  I 
First  we  shall  solve  the  sf ri njz  cigenproblem 

ti  '•  Ait  fi.  0 •  i  ■  1 

/i((li  tit  I  •  O  ill 

with  an  irregular  mesh  of  first-order  elements.  The  clement  stillness  matrix  k.  and  the  element 
mass  matrix  m.  of  the  linear  element  are 


i  n  ,, 

u  la 

,  i  •  :h 

i  i  J 

1  a  a 

When  ri  \  and  o  1,  in,  becomes  varialionally  consistent,  and  humped,  respectively.  When 
«  a,  in,  becomes  optimal  -the  accuracy  of  A  computed  with  a  uniform  mesh  of  these  elements 
increases'  from  Ot/r’i  to  Of/i  \ 

Suppose  we  are  interested  in  the  lirst  eigenvalue  A  i  n'  of  equation  1 1 1  with  the  correspond¬ 
ing  eigenfunction  n  -sin  n\.  In  the  consistent  formulation  the  energy  error'  in  the  eth  element  is 
proportional  to  hr\u"\,  where  /»,.  is  the  element  size.  Hence,  the  optimally  graded  mesh  is  here 
with  small  elements  near  the  centre  of  the  string  where  u"  is  large,  and  with  large  elements  near 
the  ends  of  the  string  where  ti"  is  small. 

To  observe  more  precisely  the  dependence  of  the  discretization  accuracy  of  A|  on  the  mesh 
grading  the  string  is  divided  into  V<-  finite  elements  symmetrically  graded  according  to 

h..  -  h,-\  e  1,2 . (.VV  -li/2  (3t 

where  e  is  changed  to  achieve  various  mesh  ratios  rj 

ri  I'a  Vc  *  ll| 

For  instance,  for  1/  2,  j  (1-38685  while  for  rj  -  ',  .* 

distributions  are  graphically  shown  in  Figure  1 . 
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The  fundamental  eigenvalue  A ,  of  the  string  is  now  computed  with  the  consistent,  lumped  and 
optimal  m,  in  equation  (2)  for  different  values  of  the  mesh  ratio  rj.  and  the  results  of  this 
computation  are  shown  in  Figure  I  As  predicted,  with  the  consistent  m,  the  optimal  mesh  ratio 
r)  is  less  than  1  and  is  seen  to  be  dose  to  one-third  It  is  interesting,  however,  that  with  our 
grading  rule  ( 3 )  the  accuracy  of  the  computed  A ,  does  not  change  much  between  17  =  3  and  q  =  J. 
A  computational  procedure  to  automatically  grade  the  mesh  is  not  likely  to  prove  profitable 
here.  We  further  note  in  Figure  1  that  with  the  lumped  mr,  grading  of  the  mesh  should  be  done  in 
the  other  direction,  with  larger  elements  near  the  centre  rather  than  near  the  ends.  Mesh  grading 
based  on  variational  arguments  will  lead  with  a  lumped  mass  matrix  to  a  loss  of  accuracy.  The 
optimal  element  mass  matrix  is  extremely  accurate  when  the  mesh  is  uniform  but  its  accuracy 
drastically  drops'  with  departure  from  mesh  regularity. 

QUADRATIC  HLFMF.NT 

Mow  do  matters  change  when  quadratic  elements  arc  used  to  discreti/e  equation  ( 1 1?  To  find  out 
we  use  the  element  matrices 


k,-l 

"  7 

8 

-8  r 
16  -8 

1 

1,  m,. 

42  r 
2  16  2 

,  m\  -  \>h 

"1 

4 

6/1 

.  1 

8  7. 

\  15 

-1  2  4  j 

1. 

for  quadratic  elements  of  size  2h.  in  equation  (5),  k,  denotes  the  element  stiffness  matrix,  inr  the 
consistent  clement  mass  matrix,  and  m'  the  corresponding  lumped  matrix. 

With  the  consistent  formulation,  the  error  in  the  eth  element  is  proportional  to  But 

here  j«'"j  -  7r'|cos  rr.ri  and  hence  the  optimal  mesh  is  with  large  consistent  elements  near  the 
centre  of  the  string  and  with  small  elements  near  the  ends. 

Because  of  their  high  order,  three  of  the  quadratic  elements  are  sufficient  for  reasonable 
accuracy  in  A They  are  symmetrically  graded,  as  shown  in  Figure  2,  and  the  mesh  ratio  tj  is 
again  increased  from  one-third  to  3.  Figure  2  verifies  the  theoretical  prediction  of  optimal  rj  >  I . 
The  gain  in  accuracy  between  q  1  and  q  3  is  greater  here  than  with  the  first -order  elements 
but  is  still  not  so  great  as  to  justify  an  expensive  computational  procedure  to  locate  the  optimal  q. 


l-si'iiic  \o  uriux  <»l  '  •  sornputcd  with  quadratic  elements 


With  the  lumped  m,.  the  grading  direction  is  once  more  opposite  to  the  consistent:  17  >  1 
causes  a  loss  in  the  computed  A  Also,  with  the  lumped  in,.,  the  error  in  A,  may  change  sign 
and  hence  the  horns  on  the  error  curve  for  this  element. 
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Problems  for  which  the  variational  principle  of  minimum  potential  energy  breaks  down,  and  are 
therefore  formally  unsuitable  for  a  finite  element  analysis,  are  shown  to  possess,  nevertheless,  a  useful 
discrete  solution. 


1.  Introduction 

Finite  elements  are  at  their  best  [1],  [2]  in  problems  that  possess  a  true  minimum  or 
maximum  variational  principle.  For  then  the  energy  convergence  of  the  finite  element  solution 
is  generally  proved,  and  the  numerical  stability  of  the  discrete  algebraic  system  set  up  with 
finite  elements  is  guaranteed  for  general  mesh  layouts.  There  are  boundary  value  problems  of 
physical  significance  and  interest  for  which  the  variational  principle  breaks  down  because  of 
boundless  energy  in  the  solution  or  because  an  analytic  solution  does  not  even  exist  for  these 
problems.  These  problems  are  theoretically  off-limits  for  the  finite  element  method  in  as  much 
as  the  theoretical  support  of  the  method  does  not  cover  them.  Energy  convergence  -  even 
pointwise  convergence  -  as  discussed  in  the  literature,  becomes  meaningless,  for  how  does  one 
measure  convergence  to  a  solution  that  does  not  exist?  Nevertheless,  conventional  finite 
element  discretization  of  such  problems  produce  discrete  solutions  that  improve,  in  some 
sense,  with  the  refinement  of  the  mesh,  and  are  therefore  entirely  useful. 

Some  simple  examples  are  solved  in  closed  form  or  numerically  in  the  following  text  to 
demonstrate  that  the  usefulness  of  finite  elements  is  retained  in  these  marginal  areas  where  the 
analytic  solution  or  the  variational  principle  totters.  Hopefully,  more  light  is  shed  thereby 
and  understanding  deepen  on  the  twin  prime  theoretical  questions  of  existence  and  con¬ 
vergence  of  the  finite  element  solution  when  viewed  from  the  discrete  point  of  view. 


2.  Circular  membrane  under  a  point  load 

Operationally  this  axisymmetrical  problem  is  formulated  as: 

-(«')'  =  0,  0<  r  <  1, 

-ru'(0)=  1,  «(1)  =  0,  (1) 
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and 


u(r)  =  log(l/r).  (2) 

The  energy  stored  in  a  unit  membrane  axisymmetrically  deflected  is 
1 

E[u]  =  ±JVrdr,  (3) 

0 

which  is  infinite  for  the  logarithmic  u  in  eq.  (2).  Let  u  now  be  any  continuous  function 
that  satisfies  the  condition  u(l)  =  0.  Formally,  the  total  potential  energy  of  the  point-loaded 
membrane  (1)  thus  deflected  is 

rr(u)  =  E[n]  -  u(0).  (4) 

To  approximately  solve  problem  (1)  with  finite  elements  by  minimizing  7t(m)  in  eq.  (4),  a 
uniform  mesh  of  N  linear  finite  elements  is  laid  upon  the  membrane.  From  E[w]  in  eq.  (3)  an 
element  stiffness  matrix  kt  is  derived  in  the  form 

*,  =  |(2e-l)[_J  “}],  e  =  1,2, . . . N,  (5) 


and  the  discrete  total  potential  energy  is  assembled  according  to  eq.  (4).  Minimization  of  ir(u) 
with  respect  to  the  nodal  variables,  under  the  restriction  that  u(l)  =  0,  produces  the  algebraic 
system  Ku~  f  for  the  nodal  unknown  vector  u,  where  the  global  stiffness  matrix  K  and  the 
load  vector  /  are 


1  -1 

f 

-1  4  -3 

-3  8  -5 

-5  12  -7 

.  /  = 

-7  16  -9 

(6) 


The  closed-form  solution  of  Ku  =  /  yields 

4.-Jj  2fTy  «'  =  1,2,...  Af,  (7) 

for  the  ith  node  and  £*+i  =  0.  At  the  center  r  =  0  (i  =  1)  of  the  membrane  where  the  point 
force  acts 


(8) 
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Even  without  a  closed-form  solution  an  assiduous  analyst  performing  numerical  computations 
with  a  succession  of  refined  finite  element  meshes  will  soon  discover  that  ut  keeps  ever 
increasing  with  N,  and  that  correspondingly  the  discrete  total  potential  energy  keeps  ever 
dropping  as  N  is  increased. 

Away  from  the  singular  point  r  =  0  pointwise  convergence  of  the  finite  element  piecewise 
linear  solution  u  to  u  =  log(l/r)  takes  place  with  the  increase  of  N.  Consider,  for  example,  the 
point  r  =  1/2  (i  =  Nl 2+  1,  N  even).  At  this  point  u(l/2)  =  log  2  and,  according  to  eq.  (7), 

“  (2)  =  2(NTT  +  7m+NT5  +  ‘"  +  N  +  N-l)- 

To  elicit  the  rate  of  convergence  of  the  series  in  eq.  (9),  it  is  helpful  to  notice  that 
1  1 

/pdp  =  l(7ydp  =  logr  (10) 

r  r 

so  that 

j  (py  dp  =  tog  2.  (11) 

1/2 

Numerical  integration  in  eq.  (11)  produces  the  series  of  eq.  (9).  Indeed,  if  |<p  <  1  is  divided 
into  N/2  equal  intervals,  then  at  each  node  thus  created  p,  =  5  +  (»  -  1)/N  (1  =  1, 2, . . .  jN  +  1). 
Within  each  interval  (p1)’  =  (N  +  2 i  -  1  )//V,  and  the  contribution  of  each  interval  to  the  integral 
in  eq.  (11)  is  2/(N  +  2 i  -  1).  Summation  readily  produces  eq.  (9).  Numerical  differentiation  and 
integration  arguments  lead  to  the  conclusion  that  this  sum  converges  0(AT2)  to  log  2.  The  same 
conclusion  could  have  been  reached  using  the  trapezoidal  rule  on  the  first  integral  of  eq.  (10). 

Fig.  1  compares  the  computed  u  with  the  exact  u  =  log(l/r)  for  N  =  10.  At  the  origin  the 
approximation  is  of  course  bad,  but  away  from  it  agreement  between  the  computed  u  and  the 
theoretical  u  rapidly  becomes  excellent. 

Once  our  ideal  analyst  has  noticed  the  singular  behavioT  of  u  near  the  origin,  he  knows  that 
he  should  change  to  a  nonuniform  mesh  of  finite  elements  for  higher  computational  efficiency. 
We  suppose  him  endowed  with  broad  knowledge  and  deep  understanding  of  the  method  of 
finite  elements  and  its  theoretical  intricacies.  Consequently,  to  decide  the  mesh  grading,  he 
sets  out  to  estimate  the  curvature  of  u.  This  he  does  by  computing  \u"\,  where 

&:  =  N2(u<-,  -  2«,  +  Ui+t),  i  =  2,3, ...  N.  (12) 

In  more  complex  situations  Q”,  may  be  computed  from  the  computer  output.  Here  eq.  (7)  leads 
to  the  formula 


A  H 

u",= 


4  N2 


(2.-1X2.-3)’ 


i  =  2, 3, ...  N, 


(13) 


indicating  the  sensibility  of  a  much  finer  mesh  nearer  r  -0  than  near  r  =  1. 
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u=ioglA 


u,  computed 


Fig.  1.  Point-loaded  circular  membrane  discretized  with  10  linear  finite  elements. 
When  the  mesh  is  nonuniform,  then  at  the  ith  node 


where  ri=0  and  rN.,  =  1,  and  in  fact  as  grading  becomes  steeper,  or  r^ilr,  ►  1,  ut  ap¬ 
proaches  a  growth  proportional  to  N.  The  rule  r,  =  «(/'  -  If,  1  =  eN",  for  instance,  yields  for 
the  specific  a  =  4  the  computed  central  deflection 


*  .15 .65  .175  ,  369,  \ 

2V  17  +  97  +  337  881  ) 


instead  of  eq.  (8). 

Energy  convergence  of  the  singular  membrane  problem  of  this  section,  excluding  the 
neighborhood  of  the  singularity,  can  be  proved  also  generally  by  the  standard  finite  element 
arguments.  Corresponding  to  IT  in  eq.  (6)  the  global  flexibility  matrix  is  given  here  by 


which  means  that  removing  the  force  from  the  origin  and  fixing  it  at  node  i  does  not  change 
the  solution  for  j  2  i,  while  Uf  remains  constant  for  /si.  The  analytical  solution  u  behaves 
the  same  way  since  -ru'  *  1  is  true  not  only  for  r  -  0  but  for  all  r.  Now  u  is  regular  inside  all 
the  finite  elements,  and  consequently  the  full  rate  of  convergence  of  the  finite  element 
approximation  to  our  original  solution,  excluding  the  singular  origin,  takes  place. 
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Table  1.  Computed  finite  element  approximation  u  to  the 
deflection  of  a  square  membrane  with  a  central  unit  force 


X 

y 

u,/V„=20 

u,  K,  =  60 

Analytical 

0.05 

0.03 

0.27531  E-2 

0.27375  E-2 

0.27357  E-2 

0.10 

0.10 

0.11015  E-l 

0.10956  E-l 

0.10949  E-l 

0.15 

0.15 

0.24840  E-l 

0.24714  E-l 

0.24699  E-l 

0.20 

0.20 

0.44456  E-t 

0.44243  E-l 

0.44217  E-l 

0.25 

0.25 

0.70505  E-l 

0.70177  E-l 

0.70138  E-l 

0.30 

0.30 

0.10451  E0 

0.10399  E0 

0.10383  E0 

0.35 

0.35 

0.14994  E0 

0.14810  E0 

0. 14890  E0 

0.40 

0.40 

0.21592  E0 

0.21333  E0 

0.21312  E0 

0.45 

0.45 

0.33650  B0 

0.32429  E0 

0.32337  E0 

0.50 

0.50 

0.72346  E0 

0.89016  E0 

X 

Such  simple  explicit  reasoning  does  not  apply  to  the  square  membrane  problem  with  a  point 
force  at  (f,  tj) 

u*x  +  uyy  +  S(£  rj)  =  0,  0<x<l,0<y<l, 

u  =  0  on  the  boundary,  (17) 

for  which  Green’s  function  [3]  is 

X,  "tSSJSL,  (l8) 

including  a  singularity  0(log  1/r)  around  (£  rj). 

To  observe  the  behavior  of  the  finite  element  solution  to  eq.  (17),  the  square  with  a  unit 
point  force  at  f  =  r)  =  \  is  discretized  with  N„  bilinear  elements,  and  the  discrete  problem  is 
solved  twice,  once  with  20  elements  per  side  (N„  =  20)  and  then  with  60  elements  per  side 
( N„  =  60).  Table  1  compares  the  computed  results  with  those  given  by  eq.  (18).  An  addition  of 
about  one  significant  digit  with  the  mesh  refinement  is  evident.  Next  the  serious  analyst  sets 
out  to  compute  u«,  uyy  and  uIy  to  better  lay  a  nonuniform  mesh  of  elements,  triangular 
perhaps,  around  the  singular  point,  but  we  are  not  yet  prepared  to  follow  him  that  far. 
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The  total  potential  energy  corresponding  to  boundary  value  problem  (19)  is 


ir(u)  =  ~j  mV  dr  -  m(0). 


where  u  is  continuous  and  u(l)  =  0.  From  ir(u)  an  element  stiffness  matrix  is  derived  in  the  form 


/t  =  3^(3e2-3e  +  l)[_J  \], 


where  N  denotes  the  number  of  equal  finite  elements  in  the  discretization.  Following  the  standard 
finite  element  procedure  produces  the  albebraic  system  Ku  =  /,  where 


K  =  — — 
3  N 


1  -1 

-1  8  -7 

-7  26  -19 

-19  56  -37 

-37  98  -61 


i)*  N%y2-y+i 


At  the  center 


“•-N%w=srrr3M5N 


since  the  sum  of  the  infinite  series  2  1/y 2  has  a  limit.  Again,  as  the  mesh  of  finite  elements  is  refined, 
the  computed  central  deflection  d,  keeps  ever  increasing,  and  with  it  the  energy  stored  in  the 
membrane. 

Away  from  the  origin,  pointwise  convergence  of  the  computed  u  to  the  exact  u  -  1/r-  1 
takes  place.  Consider  for  instance  the  point  r  =  j  (i  =  N/2+  1,  N  even)  at  which  u  =  1/(1)  -  1  = 
1.  At  this  point 


,®-"JL.t=u=v 


The  convergence  of  the  sum  in  eq.  (21)  to  u(j)  =  1  can  be  proved  by  differentiation  and 
integration  arguments,  as  has  been  done  for  the  circular  membrane  of  the  previous  section, 
from  the  equation 


r 

V.-,:  • 


5!tvpE'-' 
'  '  .r  A 


i 
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or  in  particular  from 


“(§)”  f  (Pydr- 


Division  of  the  interval  Isrsl  into  N/2  equal  segments  and  the  performance  of  the 
numerical  differentiation  of  r3  and  then  a  midpoint  numerical  integration  of  3 /(r3)'  produces 
the  sum  in  eq.  (26).  Convergence  of  u(j)  to  the  value  1  occurs  at  the  rate  of  AT2.  Fig.  2 
compares  the  computed  u  and  the  exact  u  =  1/r  -  1  when  N  =  10. 

All  the  above  discussion  assumes  a  uniform  mesh  which  is  very  inefficient  here  since 


I  “  1 

“ 7 = 18N 3  i/Mi-i)3] 


A  much  finer  mesh  near  r  =  0  is  called  for  than  near  r  =  1.  With  such  a  mesh 


fi,  =  3  r,  =  0,rN.,  =  l 

/-<  v/+i  •  r)+irj  t  rj+i 

or,  when  r;  =  e(j  -  1)"  and  cNa  =  1, 

u  =  3N“  y _ — 

■  far +r(j- it +(y-if  ’ 

and  for  large  /V 

u,  =  c(a)N“, 

where  c(l)  =  3.91,  c(2)  =  3.64,  c(3)  =  3.38,  c(4)  =  3.20  and  e(5)  =  3.10. 


u, computed 


2  8.160  9.000 

3  3.874  4.000 

4  2.295  2.333 

5  1.484  1.500 

6  0.9924  1.000 

_7 _  0.6627  0.6667 

6  0.4265  0.4286 

9  0.2489  0.2500 

10  0.1107  0.1111 


Fig.  2.  Point-loaded  spherical  membrane  discretized  with  10  linear  finite  elements. 
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4.  String  under  point  torque 

Consider  a  string  of  unit  length  with  two  equal  but  opposite  forces  P  acting  on  it  at  the 
points  j  -  e  and  \  +  «,  and  such  that  2 eP  =  M,  where  A#  is  a  constant  moment.  When  e  -* 0,  the 
deflection  of  this  string  becomes 

u(x)  =  -Mx  0  ^  x  < 

1  (33) 

u(x)  =  -M(x  -  1)  -<  x 

and  at  *  =  1  the  string  deflection  becomes  multivalued:  u({ -  0)  =  -\M  and  u(\  +  0)  =  \M.  The 
energy  needed  to  thus  deform  the  string  is  infinite. 

Our  hypothetical  analyst  ignores  this  last  unpleasant  fact  about  string  energy,  and,  being  an 
ardent  admirer  of  the  method  of  finite  elements  sets  out  to  solve  it  numerically  by  this 
technique. 

Associated  with  the  present  string  problem  is  the  total  potential  energy 
1 

ir(u)  =  u'2  djc  -  Mu'  (34) 

0 

Because  u'  is  needed  at  the  point  x-\,  cubic  C1  finite  elements  suggest  themselves  for  the 
discretization.  What  results  from  such  a  computation  is  graphically  shown  in  figs.  3  and  4  for  a 
20-element  mesh  and  M=  1.  Close  to  the  center,  where  the  torque  enters,  the  computed 
displacement  function  oscillates  severely  over  some  two  or  three  elements  (over  each  side)  as 
it  tries  to  adapt  to  the  discontinuity.  Outside  this  internal  boundary  layer  the  oscillations 

rapidly  subside,  and  u  and  u  agree  excellently.  Figs.  3  and  4  unmistakenly  suggest  a 


Hg.  3.  String  under  central  point  torque  discretized  with  20 
cubic  C1  finite  elements  (computed  and  analytical  displace¬ 
ments  shown  for  the  right  half  of  the  string). 


Fig.  4.  String  of  fig.  3  (computed  and  analytical 
slopes  shown  for  die  right  half  of  the  string). 
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nonuniform  mesh  of  finite  elements  with  larger  elements  near  the  ends  of  the  string  and  much 
smaller  elements  near  the  singular  point  x  =  |.  A  repeated  computation  with  such  a  mesh  will 
reveal  that  the  boundary  layer  is  mesh-dependent,  and,  as  the  elements  close  to  x  =  5  become 
smaller,  it  shrinks.  The  oscillations,  it  is  concluded,  are  spurious  and  are  caused  by  a 
discontinuity  of  u  at  x  =  1  with  infinite  slope  and  hence  infinite  energy. 


5.  Point  displacement  given  to  a  membrane 

Suppose  that  the  unloaded  axisymmetrical  membrane  is  given  a  unit  central  displacement. 
The  resulting  deflection  u  of  the  membrane  is  determined  by  the  solution  of 

-(ru')'  =  0,  0<r  <  1,  (35) 

u(0)=l,u(l)  =  0, 

w(r)  =  c,log(l/r)  +  c2,  (36) 

where  ct  and  c2  are  to  be  fixed  by  the  boundary  conditions  w(0)  =  1  and  u(l)  =  0.  But  it  is 
impossible  here  to  satisfy  both  of  them.  The  condition  u(l)  =  0  is  met  with  c2  =  0,  but  there  is 
no  nonzero  c ,  that  makes  u(0)=  1.  According  to  the  linear,  small  deflection  theory  of  the 
membrane  the  central  displacement  given  to  the  membrane  is  not  transmitted;  the  boundary 
value  problem  (35)  has  no  solution. 

But  what  if  one  uses  the  total  potential  energy 
1 

ir(u)  =  ij«'2rdr,  u(0)=  1,  u(l)  =  0,  (37) 

0 

to  construct  a  finite  element  solution  to  eq.  (35)?  Here  the  problem  is  not  infinite  energy.  In 
fact,  ir(u)  in  eq.  (37)  is  bounded  from  below  by  zero,  and  this  limit  can  be  approached  as 
closely  as  one  wishes  with 

u  =  1  -  ra,  a  >  0  (38) 

which  satisfies  both  end  conditions  w(0)  =  1  and  w(l)  =  0.  Introduction  of  u  of  eq.  (38)  into 
ir(«)  of  eq.  (37)  results  in 

*(»)■  (39) 

and  ir(u)~*  0  as  a  -»0. 

Standard  finite  element  approximation  of  ir(u)  in  eq.  (37)  with  a  uniform  mesh  of  linear 
elements  produces  the  linear  algebraic  system  KA=f  with 
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and 

ui  =  f,  i  =  l,2,...N-l,  «o  =  1,  nN  =  0,  (41) 

•So 

where 


i 


12/  +  1' 


At  r  =  5  (i  =  N/2,  N  even) 


, .  j  log  2 


(42) 


(43) 


and  u(j)-»0  as  N-»®. 

To  check  the  need  and  amount  of  mesh  grading,  u,  in  eq.  (41)  is  used  to  compute 
u"  =  N~2(Ui-. i  -  2u,  +  4*,);  this  yields 


A  M 

M*  = 


AT2  2 
So  4?— T’ 


i  =  1,2, . . .  N  -  1. 


(44) 


It  is  inferred  from  u"  that  a  more  efficient  finite  element  model  is  achieved  with  a  finer  mesh 
close  to  r  =  0  rather  than  close  to  r  -  1. 

Fig.  5  shows  the  computed  deflection  of  the  membrane  discretized  with  10  and  100  finite 
elements,  once  uniformly  distributed  and  then  graded  according  to  the  element  size  formula 
h,  =  «/\  It  does  not  escape  our  vigilant  analyst’s  eye  that  the  central  deflection  imparted  to  the 
membrane  does  little  to  disturb  the  rest  of  it. 


Fig.  5.  Centrally  displaced  circular  membrane.  Curve  (a)  for  N  -  100,  h,  «  j*\  (b)  for  N  -  10,  ht  «  j*\  and  curves 
(c)  and  (d)  for  a  uniform  mesh  with  N  -  100  (h  -  1/100)  and  N  -  10  (h  -  1/10),  respectively. 
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6.  Overdetermined  fixed  string  problem 

A  string  that  is  fixed  at  its  ends  is  discretized  with  cubic  C"  finite  elements.  Approximate 
solution  of  this  problem  through  the  minimization  of  the  total  potential  energy  includes  the 
requirement  that  the  finite  element  trial  function  satisfy  the  zero  displacement  provision  at 
both  ends.  Suppose  that  an  additional  constraint  is  put  on  the  finite  element  trial  function, 
namely  that  also  the  slopes  vanish  at  the  end  points.  Both  analytically  and  physically  this 
additional  boundary  condition  is  meaningless;  the  string  cannot  resist  a  torque.  The  addition 
of  the  zero  slope  conditions  merely  implies  a  slope  discontinuity  at  the  ends.  It  is  whatever 
results  from  the  correct  boundary  conditions  when  approached  from  the  interior,  but  then 
suddenly  becoming  zero  over  the  end  supports. 

To  observe  the  behavior  of  a  finite  element  solution  subject  to  the  additional  boundary 
condition  of  zero  slopes,  the  string  is  discretized  with  20  cubic  C  elements,  and  a  point  load  is 
applied  to  it  at  the  center.  At  this  point  C  continuity  is  not  enforced  on  the  trial  function  to 
allow  it  to  duplicate  the  slope  discontinuity  over  the  point  force.  At  the  end  points,  however, 
the  finite  element  trial  function  is  made  to  satisfy  both  conditions  of  zero  displacement  and 
zero  slope. 

Figs.  6  and  7  show  the  results  of  this  computation  for  half  (for  symmetry  reasons)  the  string. 
Near  the  end  point  x  -  0  the  computed  slope  u'  oscillates  wildly  inside  a  boundary  layer  of 
some  2  elements  and  then  settles  close  to  the  exact  solution  u'=  1.  The  computed  displace¬ 
ment,  on  the  other  hand,  suffers  permanent  damage  from  u'(0)  =  0,  that  is  it  is  transmitted  all 
the  way  to  the  interior  of  the  string. 

But  we  have  no  doubt  that  our  numerical  analyst,  who  is  by  now  a  seasoned  veteran  of 
many  a  difficult  computational  struggle,  will  soon  discover  that  the  boundary  layer  effect  is 
spurious  and  is  due  to  the  redundant  imposition  of  m'(0)  =  0. 


Fig.  6.  Displacements  of  fixed  string  discretized  with  10 
cubic  C1  finite  elements.  Overdetermined  boundary 
condition  of  zero  slope  at  x  -  0. 


Fig.  7.  Slopes  of  string  of  fig.  6. 
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NONLINEAR  FINITE  ELEMENT  COMPUTATION  OF 
THE  EQUILIBRIUM  AND  STABILITY  OF 
THE  CIRCULAR  PLATE + 

ISAAC  mil  l>t 

Boston  l 'niccrsttv.  Dt'ffiirftncnl  n(  Ktathcmatus.  Boston,  Massachusetts.  V.S.A. 


INTRODUCTION 

Numerical,  reduced,  integration 1  of  the  total  potential  energy  is  used  to  derive  the  element 
tangent  vectors  and  matrices  for  the  largely  deflected  circular  plate.  All  the  element  data  is 
expressed  in  terms  of  few  numerical  element  vector  and  matrices,  and  in  a  form  convenient 
for  standard  assembly  and  use  in  the  Newton- Raphson,  or  other  iterative  solution  methods. 

Actual  numerical  computations  are  carried  out  to  study  the  bending  of  the  circular  plate 
under  the  action  of  a  lateral  load  and  a  rim  thrust  that  exceeds  the  critical  value. 

A  brief  numerical  study  is  made  of  the  discretization  accuracy. 


ELEMENT  VIKTORS  AND  MATRICES 

Consider  a  unit  circular  plate  largely  bent’  under  the  action  of  a  distributed  lateral  load  /  and 
an  edge  compression  />.  1  he  total  potential  energy  of  the  deflected  plate  can  be  written  as 


rrl 


'  2  j„  (  12  (  M  *  r:  W  ')  +(‘r)  |rdr 

(«'  « >'  r  dr  -  |  fwr  dr  +  ptt(\) 


(1) 


where  «  and  «■  denote  the  inplane  and  lateral  displacements,  respectively,  and  where  (  )'  =  d/d r. 

In  the  finite  element  discretization  of  nUt,  ivl,  we  choose  to  interpolate  w  cubically  with 
the  beam  shape  functions.  Because  u  is  differentiated  only  once  in  rrlu,  it  ),  we  correspondingly 
interpolate  it  quadraticaily,  and  decide  to  numerically  integrate  rrlu,  m  )  with  two  Gauss  points, 
that  exactly  integrate  *•'"  and 

A  typical  finite  element  (see  Figure  1 1  extends  between  r  -  r ,  and  r  -  r»,  is  with  three  nodal 
points — two  end  points  and  one  central — and  is  associated  with  the  nodal  values  vector 

til  =  (Ml,  H'|,  '/riv'i.  Its,  Its,  U’l.  ’fitf  ’v  I  (2) 

where  h  =r, -r(.  To  prepare  for  the  numerical  integration  of  the  total  potential  energy  the 
typical  clement  is  mapped  to  the  standard  interval  -  1  *£  £  «T  by  r  =  Jt  1  -£)ry  +  jll  so 
that  dr  -  J/i  df.  We  denote  differentiation  with  respect  to  £  by  (  l  and  have  that  (  )'  =  2fhi  I  . 

1  Research  supported  hy  Ihc  Ollicc  of  Naval  Research  under  Contract  ONR  NOOOI4  7M‘-ll.tfi. 
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For  a  typical  interior  element  we  have  now  that 

jt.. t li.  iv )  -  |  h  'riv  ■»  ,'.r  'll  +  i/ir  '  n't  df 

t  h  '  J  ( lu'i  •  iv  r  rdf  >/t  |  /nr  df  (3) 

Interpolation  of  n  and  iv  inside  the  typical  element  is  formally  written  as 

ii  ujif  If  I  n-  -  «,!</>(£ l  (4) 

where  u,  is  the  element  nodal  values  vector  of  equation  ( 2 1  and  where  the  shape  functions 
vectors  y  and  <//  are  here 

y'tfl  ['.ftf  It.O.O,  I  f’,  'ftf  Hi.  0,0]  (5) 

and 

(/.’if) -{(0.2  3f  *•  £  \  1  f  f  ‘-f'.0,0.2»3f  f’,  1  ftf’  +  f'] 
where,  we  recall  I  •  f  *  1 . 

Next  we  substitute  a  and  n  in  equation  i4i  into  n,.lu.  h  i  in  equation  (3)  and  numerically 
integrate  it  by  sampling  the  integrand  at  the  two  C  iauss  points  f  i  s  3 '3  and  £:  =  s  3/3.  to  have 

ir,.(ttrl  V  [’,/i  V.w-;  »  ,-r,  h  '«•;  +  \hr,  u~ 

,  1 

•  h  't/ttir,  *  iv,  )  r,  \hr,f, *v,  j  (ft) 

where  the  subscript  /  -  1.2  refers  to  the  two  integration  points  £,  -  >  3/3  and  £>  =  s.3/3. 

The  values  of  u„  w,,  ti,.  iv,  and  iv,  are  computed  from  equation  (4)  as 

u,  -  ulip,  iv,  •  til •!»,  ti,  ■■■■  w.V:  iv,  -u,ih,  (7) 

with  the  numerical  clement  vectors  y„  i/-„  y„  0,  and  i h, 

y|.,  -  ,',(1  iv'3,1),  0,  4.  1  ^s. 3,  0,0) 

y  ,'.?  -=^(i2s/3.  0.  0,  *4V3,  i2s/3.  0,  0) 

i/rj  j  =  iVtO,  9  ±  4v3,  3  ±  v  3,  0,  0,  9  *  4v  3,  3±s3i  (8) 

tj,]  :  =i(0,  3,  ±v3.  0,  0.  3.  W.3) 

iji { j  =  j(0,  i  v  3,  1 1  v  3,  0,  0,  ±  v  3,  1  T  v  3) 

where  the  upper  sign  of  s/3  belongs  to  j  -  1,  (fi  ~  -  s/3/3).  and  the  lower  to  j  =  2.  (£?  =  s/3/3). 
From  equation  (ft)  we  obtain  the  element  gradient  r,  by  differentiation  with  respect  to  u, 

•.  > 

g,.  =  =  I  [  i/i  V, k>,  r  !,h  ‘r,  1  iv, tit,  +  )hr,  '«,y, 

<>u,  ,  | 

+  2/t  'r,(hu,  +  )(/iy,  *  2w^i>,)  \hr/M>,]  I9* 

since  Hu, /Hu,  -  y„  flu,/ flu,  =  y„  etc.  Or 

Rr-Y.  * a,i, +  f’l'i'i + + d&i + ‘w  •  no) 
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with 


tl. 

Fuithcr  differentiation  ol 
k.  “K-  ">• 

i  hi,  i)' n. 


\h  V, if,  />,  lh  1  r ,  'ii-,  +  4/i  V  «  Aim,  •  if •"  i 

llll 

</,  2h  r, I hit.  »■  »  :  t  r,  •  /j r ,  //, 

7r,  with  respect  t>'  u.  produces  the  element  tangent  stillness  matrix 

3,  \ii, ill,':1!  •  /’  <//,(//,'  •  (  ,<  f/r  *  i/i(e  '  i  •  </,sr,sr '  ♦  <•  ( p,*t)  |  1 1 2 1 


with 


a  ~  i/i  >r,  />,  ,’/i  V,  '  •  12/t  V,ii\  •  4  Ii  r  ti . 

«,  4 It  r, is-,  il,  2 h  c,  Wir  1 


(13) 


Now  the  element  gradient  g.  and  the  element  stillness  mntnx  k.  ate  routinely  assembled 
into  the  global  g  and  K.  the  essential  boundary  conditions  are  introduced  into  them;  the  thrust 
/),  that  appears  in  equation  tit,  is  added  to  g  at  the  entry  that  corresponds  to  nil);  and  the 
nonlinear  stillness  equation  g<»i-(>  is  iteratively  solved  for  it  with  the  N'ewton-Raphson 
method  according  to 

ir  u..  As.i'.g,,  1 1 4 1 

where  the  suhcripts  0  and  I  refei  to  the  before  and  aftei  values  ol  it  and  g. 

When  we  are  satisfied  that  the  iterative  process  in  equation  1 1 4 >  has  settled  on  an  equilibrium 
configuration  it .  we  turn  to  decide  it  stability  by  computing  the  eigenvalue  of  Kui,  i;  if  they 
are  all  positive  the  equilibrium  configuration  is  stable,  while  if  some  are  negative  the  solution 
is  unstable. 


l.ATI  RAI.  LOAD 

To  numerically  observe  the  performance  of  the  discretization  procedure  that  led  to  g,  in 
equation  110)  and  k.  in  equations  (I2f  they  are  taken  to  compute  the  deflection  of  the 
uniformly  loaded  < i  c .  /  const.),  clamped  li.e.  nil)  n  ''  1 1  Oi  circular  plate  with  an  im¬ 
movable  li.e.  it'll  i)i  edge.  Figure  I  shows  the  improvement  in  the  accuracy  in  the  central 
deflection  n  il)),  when  /  Hi.  with  the  number  of  elements  Ac.  Curve  oil  of  Figure  1  refers 
to  a  linear  interpolation  of  u,  and  <bi  to  a  quadratic.  In  both  eases  u  is  interpolated  cubieallv. 
extrapolation  to  the  limit  leads  to  the  conclusion  that  with  a  linear  u  the  relative  error  in 
(♦•iff)  is  0-46  AV  '  while  a  quadratic  u  drops  the  error  in  ivi(l|  to  (I  (108  \r  1 

Fable  I  shows  the  convergence  of  the  Newton  Raphson  method,  starting  with  the  displace¬ 
ment  vector  ii  -0,  when  f  -  Id,  ,VV  ~1.  and  the  displacements  u  and  >v  arc  interpolated 
quadratically  and  cubieallv,  respectively. 


POSTCRITIC  AL  THRUST 

For  a  unit  circular,  simply-supported  plate  of  zero  r  the  critical  rim  thrust  is  conveniently 
written4  as  p„  =  o’/ 1 2,  where  «  -  184l2(n  =  3 -99'  is  the  smallest  root  of  the  Bessel  function 
.quation  J 1  nr )  =  0 

Figure  2  describes  the  computed  behaviour  of  this  plate,  bent  by  the  combined  action  of  a 
lateral  uniformly  distributed  load  /  and  a  rim  thrust  p  that  exceeds  p,,.  When  f  >0  the  unstable 
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Table  I.  Convergence  of  the 
central  deflection  »■(())  with  the 
Newton- Raphson  method 

Cycles  n  (1)1 

1-874962 
I -.19468 1 
1  182066 
I  141467 
1142294 
1-142291 


Figure  I  Uniformly  loaded  if  |0i  tlainpeil  circular  plate  ttilli  an  immovable  idee  Computed  central  deflection 

lor  an  increased  ninnher  ol  limit-  elements  V 


a-' 


Figure  2.  Deflection-load  curve  for  a  simple  -supported  circular  plate  bent  by  a  lateral  force  f  and  an  edge  thrust  p. 
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SHORT  <  OMMl'NICATIONS 


solution  if  ()  for  p  •*  is  absent  and  a  zero  initial  guess  can  be  chosen  for  the  Ncwton- 
Raphson  solution.  The  thrust  p  is  then  increased  stepwise  with  the  last  computed  configuration 
serving  as  an  initial  guess  for  the  next  iteration.  When  f  -  0  a  non-zero  initial  guess  must  be 
used  but  the  computational  procedure  proceeds  otherwise  as  before  to  produce  the  typical 
bifurcation  curve  of  reference  5. 

At  the  bifurcation  point  p  -  p.,,  K  1  becomes  non-computablc  but  since  the  present  pro¬ 
cedure  is  global  Inon-incrementah  the  problem  of  crossing  such  a  point  does  not  arise  here. 
One  computes  the  equilibrium  configuration  of  the  plate  for  any  loading,  regardless  of  its 
history.  Fxtrapolating  from  p  >/>,,  we  find  for  f  -  0  that  near  p  -  pv„  w(0)  =  1-91(1  -  p/pc ,)I/J. 
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1.  Introduction 

The  elastica  [1,  2]  is  one  of  the  humblest  useful  elastic  systems  that  can  realistically  undergo 
very  large  displacements  with  small  strain  and  assume  multiple  stable  and  unstable  equilibrium 
configurations  under  equal  loading.  Its  total  potential  energy  is  concise,  yet  irrational  in  the 
displacements  requiring  approximate  computations  for  a  finite  element  modeling  [3-6]. 

We  find  the  elastica  a  compelling  practical  example  to  recount  the  use  of  discrete 
integration  techniques  [7]  to  derive  nonlinear  finite  elements  [8]. 


2.  Planar  elastica 

Consider  the  unit  inextensional  elastica  of  fig.  1  that  obeys  the  Bemoulli-Euler  law 

M  =  ad' ,  6'  =  dO/d s,  (1) 

which  linearly  relates  the  bending  moment  M  and  the  curvature  S’,  a  =  IE,  E  being  the 
modulus  of  elasticity  and  /  the  cross  sectional  moment  of  inertia. 

In  its  bent  state  6  =  0(s)  the  elastica  possesses  a  total  potential  energy 

ir(0)  =  JT‘  G«0'2  -  fy  +  gx)  ds  -  Py(l)  +  Qx(l)  -  M0(1) .  (2) 

Or  with  /  ■  /',  g  =  /(l)  =  J (1)  *  0,  and  since  x’  =  cos  9,  y’  »  sin  9,  x(0)  *  y(0)  =  0(0)  *  0, 

tr(9)  becomes 

*r(0)  -  £  [\a9n  +  sin  9(f  -  P)  +  cos  9(0  -  g )]  ds  -  M9(  1) .  (3) 


*Wem w>  supported  by  the  Office  of  Navel  Research  under  contract  ONR-N00014-76C-036. 
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In  this  paper  we  shall  consider  only  the  simpler  case  of 

ir(0)  =  f  (iff'2  -  P  sin  ff  +  Q  cos  0)  ds , 

Jo 


(4) 


for  which  the  admissable  ff  is  continuous  and  satisfies  the  fixed  end  condition  0(0)  =  0. 

Eq.  (3)  allows  the  expression  of  the  total  potential  energy  in  terms  of  y.  Indeed,  since 
y"  =  0’  cos  0,  eq.  (4)  may  be  written  as 

rr(y)  =  \'q  [5  +0(1-  y'T2]  ds  -  Py(l) ,  (5) 

where  y€C‘  and  y(0)  =  y'(0)  =  0. 

When  the  deflection  of  the  elastica  is  known  beforehand  to  remain  moderate  an  inter¬ 
mediate,  simpler,  theory  is  possible  based  on  the  approximations 

(l-y’2r,  =  l  +  y'*  +  y'\  (1 -y'2),ra=  1 -|y'2-sy'4  (6) 

that  change  eq.  (S)  to 

*(y)  =  1  ['  ly*U  +  y'2)  -  Qy'2(i  +  b'2)]  <*s  -  Py(  l) .  (7) 

Jo 


3.  Finite  elements 

Approximate  computation  of  the  total  potential  energy  in  eqs.  (4),  (5)  and  (7)  with  piecewise 
polynomial  interpolations  and  a  discrete  Gauss  sampling  of  the  energy  density  function  leads 
to  a  rational,  efficient  and  automatic  procedure  for  the  generation  of  high  order  nonlinear 
finite  elements. 

We  shall  first  apply  this  technique  to  ir(0)  in  eq.  (4)  which  we  intend  to  discretize  with 
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three-nodal-point,  quadratic  elements.  We  have  reason  to  believe  that  a  two  point  Gauss 
integration  of  each  element,  that  exactly  integrates  0'2  is  most  efficient  here.  A  higher 
integration  scheme  is  costlier  and  contributes  but  little  towards  improving  the  accuracy  of  the 
computed  shape.  Only  raising  the  interpolation  order  of  0  over  the  element  would  require  a 
corresponding  increase  in  the  order  of  the  integration  scheme. 

Interpolation  of  0  over  the  element  is  formally  expressed  by 

0  =  0',<t> ,  (8) 

where  0\  =  (0,,  02,  0.0  is  the  element  nodal  values  vector,  and  where 


is  the  shape  functions  vector. 

Let  the  typical  element  be  of  size  h  (not  2h )  such  that  ds  =  5/1  df,  0'  =  2 h~'0,  where 
(')  =  d/d£  Two  point  Gauss  integration  of  n(0)  for  this  element  produces  the  approximate 
element  total  potential  energy 

i 

irt  -  5/1  2  4/i ~202  -  P  sin  0/  +  Q  cos  0, ,  (10) 

/-i 

in  which  the  index  j  =  1,  2  refers  to  the  two  Gauss  points  £i  =  -V3/3  and  (2  =  V3/3 
respectively.  The  values  of  0,  and  0,  are  computed  from  eq.  (8)  as 

0/  =  0UM6).  0,  =  0'Mi) .  (li) 

We  prefer  the  briefer  notation  <f>,  =  #(£>)  and  <f>,  =  <£(£>)  and  have  from  eq.  (9)  that 
4>\  =  s(l  +  V3, 4,  1  -  V3) ,  <f> '2  =  Kl  -  V3, 4,  1  +  V3) , 

<^=!(-2V3-3,4V3,-2\/3  +  3),  <t>\  =  1(2 V3  -  3,  -4V3, 2V3, 2V3  +  3) ,  (12) 


which  we  record  once  and  for  all. 

From  nt  in  eq.  (10)  we  derive  the  element  gradient 

g«  =  =  2  Ih-'W,  -  j h(P  cos  0,  +  Q  sin  0,)<f>, 

00«  y.i 

and  the  element  tangent  stiffness  matrix  (the  element  Hessian) 

k‘  "  ^  2  2  cos  0,  -  P  sin  0,^5 , 


(13) 


(14) 


and  have  here  that 
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^!  =  T8 


<^2^2  -  Yg 


7  -8  1 

-8  16  -8 
L  1  -8  7  J 

2  +  V3  2(1 +  V3) 
2(1  +  V3)  8 

-1  2(1  -V3) 

r  2-V5  2(1 -V3) 

2(1  -  V3)  8 

-1  2(1 +  V3) 


-1  ' 
2(1  -V_3)  , 
2-V3  - 

-1  ’ 
2(1  W_3)  . 
2  +  V3  . 


(15) 


An  entirely  analogous  procedure  is  applied  to  ir(y)  in  eq.  (5),  except  that  now  y  is 
interpolated  cubically  over  the  interval  0  s  £  <  1  that  covers  an  element  of  size  h.  Once  more 
we  write  y  =  yUf>  with  the  element  nodal  values  vector  y‘.  =  (yu  yu  y2>  h)  and  with  the  element 
shape  functions  vector 

4>' « (i  -  3^ + 24U  -  2e + e,  3g2  -  n \  -e + e) ,  d6> 


from  which  <j>  and  <j>  are  computed. 

A  two  point  Gauss  integration  that  exactly  evaluates  the  integral  of  y"2  appears  to  be  most 
efficient  here  too.  These  two  points  are  in  the  interval  Osfsl  at  £i  =  1(3 - V 3),  and 
&  =  i(3  +  V 3),  and  with  equal  weights  w2  =  Consequently  n.  of  eq.  (5)  becomes 

it,  =  \h  2  \h  -  h  -2y?)-'  +  0(1  -  h  'W2 ,  (17) 

/-i 

where 

=  yi  =  y'*$i 

and 

<fc\j2  “  (— 1,  ±V3/6, 1,  + V3/6) ,  <H\,2  =  (+2V3,  - 1  +  V3,  ±2V3, 1  *  V3) ,  (18) 

where  the  upper  sign  of  V3  belongs  to  /  =  1  and  the  lower  to  /  =  2. 

From  it.  in  eq.  (17)  we  have  that 

g*  —  frtrjd y,  =  2  +  W)  >  O9) 

/-i 

with 


*>-fr-5Aft<i-*’ay?ra 


-  J*"3y/0  -  *'2y/)'‘ . 


(20) 
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a**) + mm  ) + . 


a,  =  |A-3(1  -  A-2y?r‘ ,  *>>  =  h-sM(  1  -  *'2y2r2. 

Cj  =  ih'5 y,(l  -  /r2y2)'2  +  2h-7y2,y2j(l  -  A*  2y/)-3 
-  |Qh"[(l  -  *r2y2r,/2  +  fc-2y^(l  -  At_2y>)_3/21 . 

For  the  intermediate  theory  of  eq.  (7)  we  have  in  the  same  way  that 
it,  =  J  i  h  -Jy2(l  +  h  ~2y2)  -  Oh  ' y2(l  +  3ft  '2y2) . 


gr  =  2  ^  +  bftj  - 
/-I 


a,  =  Ih-3y,(l  +  h2  y2) ,  b,  =  ^T5y2y,  -  SQh’y/O  +  ih^)-kOh^ . 


k,  =  X  +  b/(4j4>j  +  +  Cf4>i4>) . 

/-i 


a/  =  jft  ~3(1  +  ft  -2y2) ,  6,  =  h ^y,y, , 

Cj  =  Jft  -5y?  -  Sor'd  +  $/»  _2y2)  -  iOh  '3y?  • 

Notice  that 

‘  12  6  -12  6 

2  ..  _  ,  6  4-62 

I*  3E‘M/=*  =  *'  -12  -6  12  -6 

/-1  [62-64 

is  the  element  stiffness  matrix  of  the  linear  beam. 


-v  1  /A**'***  ?'&r. 


54 


/.  Fried.  Stability  and  equilibrium  of  curved  elastica 


4.  Computations 

To  locate  the  extremum  points  of  the  global  total  potential  energy  we  routinely  assemble  all 
the  element  gradients  g,  into  the  global  g,  delete  the  entries  of  g  that  correspond  to  the  fixed 
points,  add  the  tip  work  terms,  and  set  out  to  solve  the  nonlinear  stiffness  equation  g  =  0. 

If  we  choose  the  Newton-Raphson  method  for  the  iterative  solution  of  g  =  0,  then  the 
tangent  global  stiffness  matrix  K  is  also  assembled  from  all  the  element  fc,  and  the  solution 
proceeds  with 

y,  =  y„-  Ko'go,  (29) 

starting  with  some  initial  guess  that  determines  the  particular  solution  converged  upon. 

We  may  decide  that  to  update  K~l  at  each  step  is  too  expensive  and  use  instead  a  fixed  K  ' 
over  several  steps  then  update  it,  reducing  thereby  the  Newton-Raphson  method  into  one  of 
many  linear  successive  substitution  schemes.  When  and  precisely  how  to  do  all  this  is  too 
bewildering  to  contemplate  now. 

When  we  are  satisfied  that  the  iterations  have  converged  to  an  equilibrium  configuration  y« 
we  decide  its  stability  by  computing  the  (lowest)  eigenvalues  of  K(y„).  If  they  are  all  positive 
the  total  potential  energy  is  minimal  at  y.  and  this  solution  is  stable,  if,  however,  some  are 
negative  y«  lies  on  a  saddle  point  of  ir(y)  and  this  configuration  is  unstable. 

To  observe  the  actual  behavior  of  the  derived  elements  with  respect  to  the  discretization 
accuracy  and  the  performance  of  the  Newton-Raphson  method  we  undertake  to  compute  the 
deflection  of  the  elastica  given  in  eq.  (5)  and  eqs.  (18)-{21),  with  0  =  0.  Table  l  lists  the 
iterative  improvement  in  the  tip  deflection  y(l)  for  a  Newton-Raphson  computation  with 
P  =  1.5,  that  started  with  y0(s)  =  0.  Four  cycles  produce  here  a  wholly  acceptable  solution. 

When  P  is  increased  beyond  1.5  the  Newton-Raphson  method  suddenly  ceases  to  converge 
from  a  zero  initial  guess.  A  closer  initial  form  is  needed  then  to  start  the  iterative  procedure, 
or  one  may  reach  the  equilibrium  states  of  the  elastica  for  P  >  1.5  with  a  stepwise  increase  of 
P  using  the  computed  deflection  under  the  lower  load  as  an  initial  guess  for  the  next 
Newton-Raphson  iteration.  One  is  thus  presented  with  the  choice  of  small  load  increments 
with  fewer  corrections  (incremental  method  [9,  10])  or  large  load  increments  with  more 
corrections  (global  method).  The  solution  reached  with  the  incremental  procedure  depends  on 
the  chosen  loading  history,  while  the  solution  reached  with  the  global  method  depends  on  the 
initial  guess. 

Fig.  2  traces  the  computed  tip  deflection  y(l)  of  the  elastica  loaded  with  p  =  5,  Q  =  0,  as  it  is 
Table  1 

Convergence  of  the  end  deflection  y(l),  computed  with  the  Newton- 
Raphson  method,  (or  a  cantilever  beam  tip  loaded  with  a  force 
P  ■  1.5  and  discretized  with  four  cubic  elements 


Cycles 

y(D 

Cycles 

y(D 

Cycles 

y(D 

l 

0.5000000 

3 

0.4124599 

5 

0.4109928 

2 

0.4337216 

4 

0.4109994 

6 

0.4109928 
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Fig.  2.  Convergence  of  tip  deflection  y(l)  with  number.  Ne,  of  cubic  C'  finite  elements. 

improved  in  accuracy  with  the  number  of  finite  elements  Ne  used  in  the  discretization. 
Extrapolution  to  the  limit  with  the  data  in  fig.  2  discloses  that  (error  in  y(l)|  =  0.1  Ne~s  7S. 

When  the  same  elastica  problem  is  solved  with  the  element  g,  and  ke  of  eqs.  (13)  and  (14) 
we  get  for  Ne  =  3,  4,  5,  6,  7  the  corresponding  tip  angles  0  =  1.2149992,  1.2152510,  1.21531%, 
1.2153444,  1.2153549;  and  interpolation  to  the  limit  has  it  that  the  |  error  in  0(1)  |  =  0.03Ne-4. 
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Fig.  4.  Stable  equilibrium  states  of  the  elastica. 

A  seven  element  discretization  of  the  elastica  with  g,  and  kt  of  eqs.  (13)  and  (14)  is  used  to 
compute  the  stable  and  unstable  equilibrium  configurations  shown  in  figs.  3,  4  and  5.  Fig.  6 
traces  the  variation  of  Af,  the  lowest  eigenvalue  of  the  global  stiffness  matrix  K,  with  the  force 
P  for  the  equilibrium  configurations  of  figs.  3,  4,  and  5.  The  positive  branches  of  fig.  6 
correspond  to  stable  equilibrium  states,  while  the  negative  branch  corresponds  to  unstable 


**-*%>. 


Fig.  5.  Unstable  equilibrium  states  of  the  elastica. 


m 
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Fig.  6.  Lowest  eigenvalue  A*  of  the  global  stiffness  matrix  K  for  the  elastica  in:  (a)  fig.  3,  (b)  fig.  4,  6,  and  (c)  fig.  5. 


5.  Load  c  stillness  correction 

Circumstances  may  arise  in  which  the  nonlinearity  is  sufficiently  small  to  warrant  a  simpler 
successive  substitution  scheme  in  the  form  of  load  or  matrix  correction  for  the  solution  of  the 
nonlinear  stiffness  equation  instead  of  the  costlier  Newton-Raphson  scheme  with  its  involved, 
ever  reconstructed,  tangent  stiffness  matrix. 

The  intermediate  theory  of  eq.  (7)  with  0  =  0  illustrates  this.  Eqs.  (24),  (25)  and  (28)  allow 
us  to  write  the  element  stiffness  equation 

2 

ky.  =  -\h  “5  2  ffliiyifi  +  yrf>i)  (30) 

/“I 

to  be  assembled  into  the  global  stiffness  equation,  Ky  =  f(y),  K  being  the  linear,  y  free,  global 
stiffness  matrix,  and  /(y)  the  displacement  dependent  load  vector.  Successive  substitutions  is 
attempted  for  the  solution  of  Ky  =  /(y)  in  the  form  of  a  load  correction  procedure 

yi  *  K~'/(yo),  (31) 

where  K~'  need  be  formed  only  once. 

But  in  its  original  form  (31)  successive  substitutions  performs  unsatisfactorily.  To  under¬ 
stand  why  suppose  that  we  start  the  corrections  with  y0  =  0  and  compute  y  t  which  is  actually  the 
solution  to  the  completely  linearized  problem.  If  the  elastic  system  has  the  property  that  it 
becomes  stiffer  with  larger  displacements  then  y(  is  much  too  large  and  consequently  /(yt)  is 
drastically  reduced,  to  produce  in  the  next  iteration  a  too  small  y2.  Repeated  iterations 
produces  then  a  sequence  of  computed  yo,  yi, ...  that  wildly  oscillate  about  the  true  solution. 
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We  propose  to  dampen  these  oscillations  with  the  averaging  that  replaces  yi  by  (y,  +  y»)/2,  or 
change  eq.  (31)  to 

y.  =  Hyo+K'7o).  (32) 

Fig.  7  shows  the  progress  of  the  computed  tip  deflection  Y(l)  in  a  cantilever  bent  by  a  tip 
force  P  =  2,  with  successive  load  corrections,  without  (a)  and  with  (b)  averaging.  Without 
averaging  load  correction  is  useless. 


7  8  9  10 

Fig.  7.  Load  correction  solution  of  the  stiffener  equation;  (a)  without  averaging,  and  (b)  with  averaging. 


6.  Carved  elastica 

When  the  elastica  is  initially  curved  to  the  form  60  =  0o(s)  Bemoulli-Euler’s  law  of  eq.  (1) 
becomes  M  =  a{6  -  60)',  and  as  a  result  ir(6)  of  eq.  (4)  changes  into 


ir(0)  -  f  ~  9of  -  P  sin  0  +  Q  cos  d]  ds , 

JO 


or  after  integration  by  parts 


n(8)  =  f  ’  +  a$S0  -  Psin0  +  Qcos0)ds  -  a05(l)0(l) . 

Jo 


For  the  circular  arch,  in  which  03  *=  0,  ir(0)  is  reduced  to 
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Fig.  8.  Opening  of  a  circular  C-spring. 


Fig.  9.  Bending  of  a  circular  ring. 


ir(0)  =  f  ’  -  P  sin  0  +  Q  cos  0)  ds  -  a0tfl)0(l)  (35) 

Jo 

as  for  the  elastica  with  an  end  moment. 

We  use  eq.  (35),  with  or  *  1  to  compute  the  bending  of  a  circular  ring  by  two  equal  and 
opposite  forces  P,  and  the  opening  of  a  circular  C-spring  by  the  forces  O  Figs.  8  and  9  show 
the  equilibrium  states  of  the  ring  bent  under  a  sequence  of  increasing  forces.  These  figures 
compare  well  with  the  results  obtained  by  other  means  [11-14]. 
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7.  Spatial  elastica 

To  describe  the  position  and  twist  of  the  neutral  line  of  the  elastica  in  space  [1,  2]  we  need 
three  parameters  (angles)  0 ,  <p,  ift,  and  have  that 

x'  =  sin  0  cos  tft ,  y '  =  sin  0  sin  ip ,  z'  =  cos  0 ,  (36) 

where  0  =  0  corresponds  to  the  shape  of  the  free  elastica  along  the  2 -axis. 

Bernoulli-Euler's  law  becomes  here 

Mi  =  a*  i  .  M2  =  !3k2.  T  =  yr ,  (37) 

where  the  curvatures  k,.  k2  and  the  twist  r  are  given  in  terms  of  0,  <p ,  ifr  as 
k  (  =  0’  sin  <f>  -  tit'  sin  0  cos  <f> , 

k2  =  0'  cos  <f>  +  ijt'  sin  0  sin  <f> .  r  =  </>'+  ift'  cos  0  (38) 

and  the  total  potential  energy  of  the  elastica  becomes 

n(0,  <f>,  0)  =  5  f  (aK]  +  /3K:i+yr)ds-Px(l)-C>y(l)-/?2(l).  (39) 

Jo 

where  P,  Q  and  R  are  end  forces  in  the  x-,  y-  and  2-directions,  respectively. 

Finite  element  models  for  the  space  elastica  are  derived  from  ir(0 ,  <f>,  4>)  as  before. 
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Afctnct— Approximate  numerical  imefnlion  of  the  dement  total  potential  energy  with  polynomial  interpolation  of 
the  displacements  creates  hifh  order  nonlinear,  exteasMe.  cable  Inile  elements.  Successful  computations  of  static 
and  dynamic  targe  displacement  cable  problems  are  carried  out  with  the  dement. 


■rmuwcnow  y  ,p 

Gnuss  integration  of  the  element  total  potential  energy  ] 

routinely  generates  high  order  nonlinear  extensible  cable  g  /o 

inile  elements  for  any  energy  density  function  with  no  *  / 

need  for  simplifications(l).  Assembly  of  the  element  / 

matrices  and  the  Newton-Raphson  solution  of  the  global  / 

nonlinear  stiffness  equation  follows  the  liaear  inite  ele- 
went  procedure. 

High  order  elements  are  more  efficient  in  static  prob- 
Itms  and  mdispensible  in  motion [2]  problems.  Physical  ^/V\B 

discretization  with  rod  linkages  [3-6]  can  be  expected  to  L-^ S  \ _ JL 

perform  well  only  in  equilibrium  problems  of  cables  with  F%.  I.  lifted  cable 

lew  curvature. 

For  cable  dynamic  problems  with  multidirectional  ac- 

ederations  it  is  computationally  convenient  [7]  to  assume  Integration  of  cqns  (1)  yields,  with  the  boundary  con- 


independent  x  and  y  displacements  that  may  cause 
•■tension.  Near  inextensibility  is  achieved  with  a  high 
•lastic  constant. 

Large  axial  stiffness  gives  rise  to  strong  oscillations  in 
•he  computed  strain  within  each  element.  But,  as  in  the 
incompressible  elements  [S],  also  here  the  strain  and  ten¬ 
don  computed  at  the  Gauss  samplihg  points  are  accurate. 


ditions  of  eqn  (2) 


■J»+l«(l-x)+<?P 


rai—u  UUXJUU.  «<l-X>+<?  W 

_ mjgUUDSAIXOON _ 

To  lx  ideas  and  make  comparisons  we  shall  Irst  where  «-  *  Inextensibility  relates  x  and  y  to  « through 
••■potionaBy  solve  the  static  problem  of  a  balloon  lifted 

•(  the  end  of  a  long  inextensional  ideal  cable  as  in  Fig.  I.  dx  ,dr 

Anoyancy  and  wind  drag  provide  the  end  forces  P  and  Q  jjj-cosdand^-sin#  (5) 

g  the  cable  that  is  otherwise  assumed  unaffected  by  the 


5  <**•)-• 

£<#  caefl+jr-l 


•««<l 


p  denotes  the  mass  par  unit  length,  and  whan 
£■# U)  k  the  tension.  At  the  elevated  end  of  the  cable 


ten#*  AQ 

_  ?-r+<r- 

Aintd  by  d»  Office  ef  Novel  teieareh 

®5^moi4-iic-dC9*. 


To  prepare  for  a 


sre  naad  the  total 


•<AP)--nJ^  xdi  +  PXD-Qbrfl) 


1U 


)K 

of  the  incxtensibk  cable,  which  with  d(xy)  -  z  dr  + 1  dx 
and  with  eqn  (S)  can  be  written  in  terms  of  §  only: 

•(#)■•£  (r-l)eos#d*-j£  (P  m  Q + Q cos#) dx 


An  inextensional  cable  element,  with  f  nodal  values 
only,  can  be  readily  derived  from  off)  in  eqn  (9). 
However,  with  an  eye  on  dynamic  problems  and  for 
greater  generality  we  prefer  to  allow  axial  stretching  and 
dm  corresponding  additional  elastic  energy. 

In  terms  of  x  and  y  the  tensile  strain  is 

e«(x*+yT*-l  (10) 


with  /■  1,2  referring  to  the  two  Gauss  points  (,  -  -y/3/j 
and  ft  »  v/3/3.  The  values  of  a,  y,  i,  and  j  at  the  Gatm 
points,  namely  x,,  *  ^  and  A  needed  in  eqn  (17)  me 
computed  from  eqns  (UH13)  as  x. «  u/ey,  y,mx,T± 
4,-a/*andA-a/*  where*- offend* ■ 
etc.  Actually  4 

P.a-J  (1±V3.0,4.0.UV3.0) 
0.4-g(M±V3,0.4,0,l*V3) 
Pu“g(»2V3-3.0.±4V3.0,?2V3+3.0)  (0} 

♦.4-7(0,  rtV3-3,0.±4V3,«,«V3+3) 


where  (fm  d/dr,  and  for  Hook  material  the  extra  term 

ofx.yj-lcj^dr  (11) 

is  added  to  the  total  potential  energy  in  eqn  (I).  Addition 
of  the  elastic  energy  in  eqn  (11)  permits  an  axial  elon¬ 
gation  of  about  1/e. 


m  which  the  upper  sign  of  V3  belongs  to  /« 1  and  the 
lower  to  /  «  2. 

The  element  gradient  g,  and  the  element-  tavern 
stiffness  matrix  k.  are  derived  from  the  element  lot* 
potential  energy  «r,  in  eqn  (17)  through  repealed 
differentiation  with  respect  to  the  element  nodal  vrieer 
vector  a,: 


a*  TO 


We  now  describe  the  development  of  a  three-nodal 
point  cable  element  with  a  quadratic-quadratic  inter¬ 
polation  of  x  and  y  over  it,  and  a  two  Gann  point 
integration  of  the  element  total  potential  energy.  Inside 
each  element  the  coordinates  of  the  displaced  cable,  x 
and  y  are  interpolated  from  the  corresponding  nodal 
values  by 


x  »  a/p  and  y  *  a.rd  (12) 

where  the  nodal  values  vector 

•UT-(x„yl,Xj,yJ,*»,ys)  (13) 

and  where  die  Ingrange  interpolation  (shape)  functions 
are 


br-|U(*-l).0.2(l-<,).0. «(+!),  0)  (14) 

♦T«jw>,M-l),0.2(l-f>.0.«f+l>l  (IS) 


lor -1*|*1. 

With  *-(,+#,  dr-hd£  x'-fc-'i,  and 
where  0>d/d£  the  total  potential  energy  of  a  typical 
element  becomes 


*  —  «*  f  \  xd*+{e*  J’  th-'fi**  f)'n- l)d( 

(14) 

A  two  point  Gauss  integration  of  v,  in  eqn  (14)  requires 
We  evil— don  of  the  energy  density  function  at  (■ 
*V3/3  and  its  summation  with  eq— 1  unit  weights.  Hence 
o.  of  eqn  (14)  becomes  ippreximataly 

(17) 


so  that 


where 


fc-h^-W+cG-*;^  | 

$-*?+*?  — d^-xjpj+y*,*’,  01)7 

K  -  ch  i  (l  -  «7,aXp + Itf)-  «7*awT. 

(28  * 


In  the  fnite  element  computodoo  of  the  cable's  < 
fcrium  an  initial  configuration  guess  vector  *  is  made,g,  j 
and  k,  arc  computed  with  this  guess  and  are  assent—  ] 
into  die  global  $,  and  K*.  Essential  boundary  i 
(here  x(0)  -  y(0)  -  0)  art  imposed  upon  u  and  A*  as  al 
the  linear  case  and  an  improved  conhguration  a.  b  j 
1  with  the  Newtoo-Raphson  method  i 


to  a, -*-*5'*  Wj 

to  aw  Stability  of  the  computed  eqdl 


defotim—m  of  XL. 

IS— re  2  shows  the  comouled  —hie  coalnunbeo  ef 
Pi«Tl  for  <?-0  and "V- 081,  OuSTTw/lA 
0i0,...,0.90,  for  a  seven  element  dhtniireri—  erf 
c  - 100.  When  P  is  smal  the  cable  tip  behavior  is  u— 0 
siagniar  since  it  must  rotate  almost  localy  to  O-nOh 
order  to  tenge  misty  most  the  vortical  forte.  An  on*; 
dp  curvature  k  created  foereby  (winch  m#tfe 


•I  the  dp) 


Or  numerically  integrated 


».--*“**  jE  Q  y’+«a)  05) 

Consequently  an  additional  element  gradient  and  ele¬ 
ment  stiffness  matrix 

t, (A + e)6  (26) 


a.*iii4ffc"f'jr’gb 

d,-dt-}fir,(*4ff.>  (2ft 

for  a  time  interval  r. 

A  lumped  element  mass  matrix  [2]  for  the  nodal  values 
arrangement  in  eqn  (13) 


1  -f*1*  J)  m,t 


need  be  taken  into  account 
Figure  5  shows  a  series  of  computed,  unstable,  equili¬ 
brium  configurations  of  the  string  (p  ■  1)  for  *  *  0.02  and 
«•*»  15i 4 2,  i - 1, 2, 3, .... 20;  discretization  being  done 
with  seven  elements. 

CASLt  DVtUMKS 

With  «  denoting  the  cable’s  global  vector  of  nodal 
unknowns  its  equation  of  motion  is  written  as 

gfn)4Mg-0  (28) 

in  which  g  is  the  global  total  potential  energy  gradient, 
and  tH  the  global  mass  matrix.  A  step  by -step  integration 
of  eqn  (28)  starting  with  the  initial  position  a.  and 
velocity  li.  is  performed  with  Newmark’s  method  [2], 


Ne*7 
e <400 
•*25j*2 
ftZr-.ZO 


••002 


Fig.  j.  (Instable  coalgantioes  of  inrntTriag  cable. 


iu.-fdiag.Cl.  1.4.4, 1, 11  oi, 

is  used  to  assemble  M 

As  an  example  we  shall  compute  the  free  fall  of  a 
chain  originally  fixed  between  two  horizontal  points  aatf 
then  released  at  one  end.  At  rest  the  shape  of  the  hai^ 
chain  is  parametrically  given  by 


,L  t^M} 


and  j 

Our  particular  case  is  for  uf(2?)- 1,  or  x(l/2)«O201 
and  yd)-0J81. 

Figure  6  traces  the  position  of  the  chain  at  each  1/lOof 
a  second  for  computations  carried  on  with  r-  1/10*0, 
N. » 15,  and  c*  100.  Figure  7  shows  the  same  falag 
chain  traced  from  a  video  camera  Urn  of  an  actual  chsia 
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F*.  7.  Experimental  counterpart  to  Fig  *. 


tea  taken  at  not  equally  spaced  time  intervals.  It  is 
Markable  that  the  computation  agrees  with  reality  even 
thout  such  fine  details  as  the  curling  up  of  the  free  end 
takes  place  while  the  chain  flings  back  across  the 
Mkallme. 
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FINITE  ELEMENT  COMPUTATION  OF  LARGE 
RUBBER  MEMBRANE  DEFORMATIONS 
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Boston  University,  Department  of  Mathematics,  Boston,  Massachusetts,  U.S.A. 


INTRODUCTION 

The  stretched  and  inflated  rubber  membrane  problem1,2  abounds  in  realistic  examples  of 
elastic  systems  that  may  undergo  extremely  large  deformations  far  beyond  anything  a  linear 
theory  can  handle.  Geometric  and  material  nonlinearities  sufficiently  complicate  the  situation 
for  the  analytic  answer  to  its  equilibrium  question  to  become  impossible  even  for  membranes 
of  the  simplest  geometry. 

Rubber  is  conveniently  characterized  by  an  energy  density  function  and  it  comes  natural 
to  apply  variational  methods, 2,4  including  finite  elements,510  to  the  approximate  computation 
of  its  deformation. 

Piecewise  polynomial  approximation  of  the  displacements,  coupled  with  an  element-by¬ 
element  discrete  integration  of  the  total  potential  energy,  promises  to  be  the  most  general 
and  efficient  technique  for  the  solution  of  the  rubber  membrane  problem.  Being  a  finite 
element  method  this  solution  technique  is  highly  programmable  and  includes  accuracy  and 
efficiency  controls  through  high  order  elements  and  a  finer  mesh. 

In  this  paper  we  derive  in  detail  the  element  gradient  and  element  tangent  stiffness  matrix 
for  the  axisymmetric,  Mooney,  rubber  membrane,  including  a  quadratic-quadratic  interpola¬ 
tion  of  the  displacements  and  a  two-point  Gauss  integration  of  the  element  total  potential 
energy.  Such  a  discretization  procedure  is,  evidently,  indifferent  to  the  complexity  of  the 
energy  density  function  and  may  be  extensively  applied  to  other11"13  than  Mooney  materials. 

We  employ  the  quadratic  element  to  compute  the  inflated  and  stretched  shapes  of  the  disc, 
the  torus,  and  the  tube,  for  which  other  comparative  computational  and  experimental  results 
are  plentiful.14"2*  These  numerical  examples  are  made  to  check  the  correctness  of  the  formula¬ 
tion,  to  exhibit  the  versatility  of  the  finite  element  technique,  and  to  study  the  accuracy  of 
the  element.  The  convergence  of  the  Newton-Raphson  method  near  a  critical  point  is 
scrutinized. 


AXISYMMETRIC  MEMBRANE 

With  reference  to  Figure  1,  let  the  generating  curve  of  the  undeformed  membrane  be  described 
in  the  (r,  z)  plane  through  r  =  r(r)  and  z  -  z(i),  where  s  denotes  arc  length.  Under  the  action 
of  applied  forces  and  prescribed  displacements  the  point  (r,  z)  moves  to  the  deformed  location 
(x,y).  An  arc  element,  originally  dr  is  stretched  thereby  to  dr*,  and  the  membrane  thickness 
t  shrinks  to  r*. 


t  Profewor.  Rewardi  supported  by  tfcsQgae  o«  Nsvtl  Keswroa  CmiuacrOWR-N00014^eO0te<. 
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I.  FRIED 


Figure  1.  Undeformed  and  deformed  element  of  arc  on  the  generating  curve  of  an  axisymmetrical  membrane 


The  energy  density  function  of  the  rubber  membrane  is  ultimately  expressed  in  terms  of 
the  three  principal  stretch  ratios 


A  i  — 


ds* 


A2  — : 


2  irx 


A3 - 


d s '  2m'  '  t 

and  if  the  deformation  is  volume  preserving 

A1A2A3 —  1 

Since  ds*  =  (dx2 +dy2)l/\  equations  (1)  and  (2)  become 


A,  =  (x'2  +  y'2)I/2, 


A2  —  . 

T 


A3  = 


A1A2 


(1) 

(2) 

(3) 


where  (  )'  =  d/dx. 

Here  we  follow  the  common,  rather  realistic,  assumption  of  an  incompressible  Mooney 
membrane,  inflated  under  the  pressure  p,  for  which  the  total  potential  energy  is  of  the  form 

tr(x,y)  =  2tr^j{ot(/i-3)+a(/2-3)]rds  +  i^|oxVdsj  (4) 


where  p  and  a  are  material  constants  (when  a  =  0  the  material  is  modernistically  named 
neo-Hookean),  and  where  /|  and  f2  are  the  strain  in  variants 


/i=A2  +  a|  +  A3  1 

/*=a?a!+aIa!+a^ 


(5) 


Henceforth  we  shall  replace,  for  typographical  brevity,  p/pt  by  p,  and  we  shall  take  2-trpt  *  1. 


FINITE  ELEMENTS 

A  typical  three-nodal-point  element  is  shown  is  Figure  1.  Inside  each  such  dement  x  and  y 
are  interpolated  quadratically  by 

x  =  u  It,  y  =  nT*  (6) 
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where  the  element  nodal  values  vector  wj  *=  (xt,  y\,  x2,  ya,  xj,  yj),  and  where  the  shape  function 
vectors  <t>  and  #  are  here 


*T =&(*  - 1),  o,  i -e,  o,  Je(f + 1).  o] 

(7) 

*T = [o,  k(f  - 1),  o,  i  -  o,  if  if +i)] 

in  which  -1  =£ 1. 

From  s=s2  +  hf  it  follows  that  ds  =  h  df,  or  x'  =  h~lx,  y'  =  h~ly,  where  ()  =  d/df.  Two- 
point  Gauss  integration  over  each  element  of  the  approximate  ir(x,  y)  is  sufficient  to  remove 
from  the  element  all  spurious  mechanisms  and  artificial  instabilities  while  assuring  high 
computational  efficiency.  In  the  interval  -1  <s  (  « 1  the  two  Gauss  points  are  at  fi  =  -V(3)/3 
and  f2  =  V(3)/3,  with  the  equal  weights  wj  =  w2  =  1. 

The  element  total  potential  energy  irt  is  expressed  in  terms  of  seven  element  integrals 

ire  =  J\+J2+J2  +  <x(J4+Js+J(,)  +  2pJy  (8) 

numerically  integrated  as 

/.=  [  A ds  =  h~l  f  (x2  +  y2)r df  =  /i_1  £  r,(x?+y?) 
i,  J-i  /-i 

J2  =  f  \22rds  =  h  f  r  lx2df  =  h  £  rj'x2 
i,  J-i  i- 1 

/3=  f  Xl2\-22rds=h3  f'  r3(x2  +  y2)~lx~2  df  =  h3  £  +yj)~l 

Jr  J-i  i-i 

f  A^lrds  =  h~l  f  r"lxz(il +  yl)df  £  ri'xf(xf+yf)  (9) 

l,  J-i  /-i 

7S=  f  \?rds  =  h3  f  r(i2  +  y2r,<4  =  h3  £  r,(i,2  +  y?)"1 
J*  J-i  i- 1 

/6  «  f  AJ2r  di  =  h  f  x'V d£  =  h  £  r3txj2 

J*  J-l  /-I 

A  »  f  * V  dj  “  [  x2j)df  =  I  x2y2 
J,  J-i  /-I 

where  the  subscript  /  *  1,2  refers  to  the  two  Gauss  points  f2  =  ~V(3)/3  and  £>  =*  7(3)/3.  Also, 
in  equation  (9)  rt  <*  rfs/)  arid 

xi  =  uj<t>„  y,  =  uJiA/,  x,  =  i<77>,  y/  =  nJ7/  (10) 

after  setting  ^  =  ^(^).  From  equation  (7)  we  get  the  numerical  vectors 

-  J(1  ±  7(3),  0, 4, 0, 1 *  7(3),  0) 

*,.2  -  l(*27(3)  -  3, 0,  ±47(3),  0,  *27(3) +3,0)  (11) 

+a,2  «  *27(3) -3, 0,  ±47(3),  0,  *27(3)+ 3) 

where  the  upper  sign  of  7(3)  belongs  to  Gauss  point  1  and  the  lower  sign  of  7(3)  belongs  to 
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H>.  4*"  :  . 

«-'••••••  t 


Differentiation  of  7r„  given  by  equations  (8)  and  (9),  with  respect  to  ue  produces  the  element 
gradient 

ri'tr  2  t  . 

g<=-rs=  I  aA+b/ti  +  crf,  (12) 

due  i-i 

with 

%  =  2r/xJ(l  —  A  i^A  if  )(1  +  oA  |/) 

b,  =2riy',(\  -Ai?Al,2)(l  +«Al/)+!pjr?  (13) 

Cy =  2/»A2/(l  -'Ai/2A2/4)(1  +<*A  iy)  +  hpxyy} 

Further  differentiation  of  g,  with  respect  to  u,  produces  the  element  tangent  stiffness  matrix 

^  =  ^  =  ~r=  I  +  b,Mj  +  C/Ml 

dUg  OUg  fm\ 


+  +  Ml )  +  )  +/,(<£, +  Ml )  (14) 


where 


a,  =  2h  ‘r,(  1  - -'A^-/'  )(l  +  «A |,) 

bf  =  2h  1  — ^  +  “A 

Cj  —  2hfj 1  (1  +  3A  i,2A  2,4  )(1  +  aA?/)  +  hpy ) 
dj  =  4xjA2/(a  A-  A  i/A  2* ) 
ei  ~ 4y/A2,(a  +Ai,4A2 *)+pxt 
fi  —  Sh  'ryXyyjA  1  *A 2/* (1  +<*A2/) 

The  element  gradient  g,  and  the  element  tangent  stiffness  matrix  K,  are  linear  combinations 
of  the  numerical  vectors  <j>h  <L  and  and  the  numerical  matrices  6f6j,  Ml,  4>id>J,  <bi<bj  + 
Ml  +  Ml  and  <btii>l  +  Mi-  To  compute  the  displacement  dependent  coefficients 
ah  bh  ch  dh  eh  fi  of  these  combinations,  u,  is  picked  out  from  the  global  displacement  vector  u 
and  is  introduced  into  equation  (10)  to  yield,  again  with  the  aid  of  the  numerical  vectors  <)>h 
<j>i  and  \jih  the  values  of  xy,  yh  Xy  and  yy.  These  values  are  used  to  compute  the  stretch  ratios 
Ai/  =  X/2  +  yJ2  and  kh-xf/rj,  where  xj  =  h~xxt  and  y|  ~/r,y/»  with  which  the  coefficients 
in  equations  (13)  and  (15)  are  finally  computed. 

Once  ge  and  kt  are  formed  the  finite  element  assembly  procedure  follows  for  the  nonlinear 
case  precisely  as  in  the  linear:  an  initial  global  displacement  u0  is  made,  all  g,  and  ke  are 
computed  for  it  and  routinely  assembled  into  the  global  g  and  K,  the  essential  boundary 
conditions  are  introduced  into  g  and  K,  and  no  is  improved  into  ut  with  the  Newton-Raphson 
method 

ut~Ho-K~\u0)g{Uo)  (16) 

until  convergence  to  u«. 

The  lower  eigenvalue  spectrum  of  K{um)  indicates  the  stability  of  the  onmputed  aohitiou. 
AH  positive  eigenvalues  mean  that  u,x,  is  at  a  minimum  potat  of  the  total  potential  energy  and 
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is  therefore  stable.  Some  negative  eigenvalues  mean  a  saddle  point  of  ir(x,  y)  and  unstable 
equilibrium. 


COMPUTATIONS 

Finite  elements  prove  themselves  best  through  hard  work.  We  put  our  element  to  the  work 
of  computing  the  deformations  of  a  disc,  a  toms  and  a  tube. 

First  the  disc.  It  is  originally  described  by  r(r)  =  s,  z(s)  =  0, «  1,  its  edge  is  stretched 
to  x  =  1-1,  y  =  (0),  and  a  pressure  p  is  applied  to  its  face.  The  purpose  of  this  initial  stretching 
is  to  endow  the  membrane  with  a  linear  solution  and  consequently  an  easier  application  of 
the  Newton-Raphson  method  with  a  zero  initial  guess. 

At  the  outset  we  assess  the  accuracy  that  our  element  can  provide  in  order  for  the  future 
computations  to  be  correct  but  not  overly  expensive.  We  compute  the  polar  rise  y(0)  of  the 
disc  ( a  =01)  caused  by  a  pressure  p  =  5.  The  disc  is  substantially  deformed  under  this  value 
of  p,  and  for  a  uniform  mesh  of  Ne  finite  elements  we  have  that  corresponding  to  Ne  = 
1,2,3,4,10,  y(0)  =  1-2546,  1-4278,  1-4299,  1-4304,  1-4304.  Five  finite  elements  assure 
reasonable  accuracy  in  the  displacements. 

Our  next  concern  is  with  the  performance  of  the  Newton-Raphson  method.  To  form  an 
idea  as  to  how  this  method  works  we  compute  y(0)  that  results  from  p  =  5,  with  an  initial 
zero  guess  and  a  ten-element  discretization.  Newton-Raphson’s  method  successively  computes 
y(0)  =  2-31875, 1-44782,  1-44302,  1-43059,  1-43040,  1-43040;  and  four  cycles  are  sufficient 
for  a  six-digit  accuracy.  These  computations  are  carried  out  in  a  single  precision  with  some 
six  significant  digits.  Double  precision  could  have  saved  us  one  cycle. 

Close  to  a  critical  point  where  K  is  singular  the  Newton-Raphson  method  slows  down,  as 
we  shall  see  soon. 

Figure  2  shows  the  inflated  disc  for  a  pressure  that  increases  at  a  step  of  one  between  p  =  1 
and  p-  7.  Corresponding  to  these  pressures  are  the  polar  stretch  ratios  A0  =  1-144,  1-233, 
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Loss  of  stability  in  a  stretched  and  inflated  right  tube  (a  =  0)  through  bulging  is  shown  in 
Figure  4.  Simulation  of  the  bulging  is  achieved  by  pulling  out  the.  central  circumference  of 
the  tube  and  holding  it  fixed.  A  pressure  is  then  introduced  into  the  tube  that  is  left  to  inflate 
until  the  slope  at  the  central  grip  point  becomes  parallel  to  the  y-axis. 
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1.  !  N’TNOIHCTI  ON 

Approx  I  mate  Cans:,  quadrature  of  t!u'  total  potential  oncve.yi  2  I 
is  show i ns  ijroat  promise  to  become  a  universal  means  lor  the  set 
up  o(  nonlinear  (inn.  cleiv.mls.  As  in  t  lie  linear  case  also 
here  all  the  numorie.il  Iv  inlee  rated  linite  elements  ar<-  expres¬ 
sed  in  terms  of  l<w  numerical  vectors  and  matrices  and  in  a  form 
enliven  lent  for  standard  assemble  and  use  in  the  Newton-Naphson 
met  hod . 

Detailed  derivation  and  a-  t  ua  I  •  -input  at  inn  is  included  in 
this  paper  tor  the  clement  ,'radient  and  stiffness  matrices  u! 
the  1  are, el  v  deformed  bean,  rine,  circular  plain  and  rubber  mem¬ 
brane  . 

2.  h i : a: !  a:. h  k  i  . 

Cur  start  in;;  point,  for  the  unit  clast  ica  !X  I,  shown  in  i'ii».  1, 
is  its  total  potent  ini  enere.v 

i 

K  >  (.'ll  F  in  t  i>  cos  M  (  11  ( 2 .  I  ) 

n 

for  which  the  admissible  is  emit  i minus  and  satis)  ics  Lhe  1  ixed 
end  condition  •(<))  =  <). 

Initial  eurvat  ire  in  the  lorn  =»  ( s)  alters  the  total  po¬ 

tential  energy  ul  (2.1)  into 

»(■■>  =  I  l5KI(-  '-■’)•  -  I'sin  +  ('  cos  Ids  (2.2) 

’  0 

where  the  end  moment  M  is  assumed  absent.  Alter  integration  lie 
parts 
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!•'!(!.  1.  Tip  load  oil  clast  ica. 


n(»)  = 


•  2 

+  !:i  " 

I’sin-  i  ■•os 


Ids-CI  '<  1  )  <  n 


{ fl 


and  for  the  circular  rim-,  "  -  0  so  that 

»(•■)  =  I  1  f 4l  1  *-!’« in"  f  n  cos  )ds  -  KI  '  fl)  (II  (?..'<) 
'  0 

I 

.is  for  the  straight  Hast  ica  with  an  end  moment  M*Ff*  M). 


Through  the  use  ot  '  cii.i  ■  and  v' 
tal  potential  energy  in  terms  of  v  as 
'  1 


we  ma v  write  the  i  < 


tr(v) 


I  Iki 


'  (1 


I  ■ 


(  (.1(1  -v'  )  ■'  |<!s  -  I'vf  |  ) 


where  t  tic  admissible  y  is  C  and  satisfies  the  fixed  end  mo¬ 
di  l  i  oils  y  ( fl)  y '  (  0)  -  (). 


We  propose  a  finite  element  discretization  ol  it  (-  )  in  (J.l) 
with  a  quadrat  ic.  interpolat  ion  of  over  each  element  and  a  two 
point  Cans::  quadrature  ol  the  total  potential  enerr.v.  Ttiis  min 
imal  integration  scheme  assures  t  lie  numerical  stability  of  t  In¬ 
finite  element  method  and  is  snl fir  iently  accurate.  The  resul- 
tinp  quadraL  ic  element  is  precise,  efficient  and  easily  propram 
mab I e. 


interpolat  ion  ol  •  oyer  tin-  l  hree-utoda  I  -po  i  n  t  element  is  com¬ 
part  I  v  expressed  bv  "  1 1 , .  e  where  u,.  I‘ 1  ,  '  i  •  the  cle¬ 

ment  nodal  values  veclot  ,  aud 
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After  r.t.  and  VK.  are  prt  pi  red  t  he  •  inite  element  a  ssemblv  and 
solution  prmvilurr  to!  1  •  *  w  ■ .  toi  t  he  i***ti !  in*-.n  imm'  ;*n<  i.il-  as 
tut  tin*  linear:  An  in  i  t  i  1 1  •  uM  is  made  lor  tin  "  I  oha  1 

slopes  vector  u,  a  1  1  and  Vt.  are  eor.-puted  I  row  it  and  routine¬ 
ly  assembled  into  the  p.lohal  >•.  and  r\  t  iie  e.seiltial  boundary 
t  audit  ions  am!  boundary  woi  !•  teim  .  are  itn  htded  in  .•  ami  i  :  and 
u0  is  improved  into  ii]  v.ith  t  hi-  Newt  on- haphsou  method 
it,  =  ii  -  K  .•  nnt  il  sonvery.i  n*  e . 

To  discretize  »fv)  ill  (  .  r>  »  We  piop.ea  a  p  i  e«  •  -V  i  *.e  cubic,  (  , 
interpolation  oj  v  witii  I  wa  *.  ms  .  point  .(uadrature  ui  the  rle- 
!!ient  tot  ill  potential  enerpy. 

Now  y  *  uc.  ;  ,  when-  u,,  • v( ,v(  ,v. ,  and 


f,--*  -?4V\,  V* 


i  on?) 


and  I  lie  !  we  dan:,:,  pninl:  if"  at  |  <  •-■  ((/'•  and  ,  i  ,1  ')/'■ 

with  erjna !  weights  =  w  f  -  .  'ei  »rpif*ntl"  t  ae  appro:*,  imal  i 

element,  total  pot  ent  ial  emr.".  1.'.  a)  hei  rr. 
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tbi'  Newt  on-l;.iphson  method  sue erm  i  vcl  y  computes  u  tip  dot  led  i 
y(l  )--0.  5000000,  0.4  ('5721  8,  ().41245''9,  (I.- 'll  O'KKt’i ,  0.410«928, 

0 .  4  1 09928 ;  h.iviii)’  started  with  .1  zero  initial  tar,. 


When  2  is  increased  above  1.3  the  hewton-ttaphson  method  sud¬ 
denly  fails  Lo  convert’, e  from  a  zero  initial  cucss.  A  hotter 
■■tart  inp,  shape  is  needed  then  lor  the  iterative  solution,  or 
the  deflection  under  higher  loads  can  he  reached  stepwise  with 
the  computed  solution  at  tin’  end  of  the  previous  step  servim; 
as  an  init  ial  ipiess  for  the  next  iteration  with  the  higher  load, 
one  is  thus  confronted  with  the  choice  ol  small  loud  increments 
witli  fewer  iterative  correct,  ions-an  incremental  method,  or 
large  load  increments  with  more  correct  ions-a  global  method. 

In  l  lie  presence  ol  mult  iplc  solut  ions  to  tin-  stillness  eipiat  ion 
g(ti)=*0,  the  load  history  of  the  incremental  method  is  what  de¬ 
termines  which  one  of  them  will  he  discovered,  while  the  solu¬ 
tion  reached  with  the  global  method  is  determined  by  the  ini¬ 
tial  guess  Uj,. 

To  disclose  the  diserct  izal  ion  accuracy  ol  the  element  in 
( 2 . 1  5)  -  ( 2 . 1 8)  a  varying  number,  Nc,  of  elements  are  employed  in 
the  computation  of  cantilever  deflected  by  I*  -  3.  Tor 
Ne  *  1, 2,1,..., 9  we,  respectively,  compute  y(l)  =  0.7hb9329, 
0.7183933,  0 . 7 1 4 1 3 1 A ,  0.7 1  HOI  74 ,  0.71)9174,  0. 7138022,  0. 7138340, 
0.7138185  and  extrapolation  to  the  limit  provides  the  estimate 
error  in  y  ( I  )  ,  =  O.INe"*’ 

Solving  tlie  same  cantilever  problem  with  the  element  given  in 
(2.9)  and  (2.10)  we  get  Iur  Ne  =  3,4, 3,  (>,  7  the  corresponding  lip 
Slopes  (I)  -  1.2149992,  1.2152510,  1.2153190,  1.2153444, 

1  .2153549;  and  with  extrapolation  to  the  limit  we  reach  the  es¬ 
timate  ierror  In  i(l)  1  =  0.03No"  . 

A  seven  element  discretization  of  the  elastic;!  with  the  cle¬ 
ment  data  in  (2.9)  and  (2.10)  computes  the  stable  and  unstable 
equilibrium  conf  igurat  ions  14  I  shown  in  Figs.  2,3  and  4.  Figure 


I.  FRIED 


FICI.  f> .  Squecz  ing  of  a  circular  ring. 


3.  CIKCIJI.AK  P1.ATK 


A  tin  it  c  i  nil  I  a  r  plate  ( 


1 1  i  i\ 1 1 1  v  I. 


I  1  louder 


tion  of  an  ax  isvmmet  r  icu  1 1  v  distributed  lateral  load  I 


uniform  edge  compression  p. 
plate  reads 


The  total  potential  energy  of  this 


ir(u,w)  = 


fii ' +iw'")-r.lr  - 


twrdr  f  pu(l) 


where  u  and  w  denote  the  inplane  and  lateral  displacements,  res¬ 
pectively,  and  where  (  )'  =  d/dr. 

We  propose  to  discretize  f(u,w)  with  a  d'  piecewise  cubic  in¬ 
terpolation  of  w,  a  piecewise  quadratic  interpolation  of  u,  and 
a  two  dauss  point  fill  cp.ral  ion  of  I  lie  total  potential  energy  over 
each  element.  A  linear  interpolation  scheme  lor  ii  is  noticed  to 
produce  a  decidedly  inferior  element. 

Typically  un  element  extends  between  r»r,  and  r«r„,  has 
three  nodal  pyints  and  is  .associated  with  I  lie  element  nodal  va¬ 
lues  vector  u(,  ~  !ii|,W|,w  u,,  ,u  ( ,  w  ,  w.^  .  I  ul  erpo  I  at  iyu  ol  u 

and  w  Inside  I  lie  element  ift  Tormallv  written  as  u  u(,.'f  )  and 
w  ii,1..  (  )  with 


-J 
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J>>  ^',('.+D ,o,o! 

»  i  i  *>  i  ^  7 

2  “  ^{0,2-3c.+  ,  1 ,  0,0, 2+3-—  ,-!->•  +  ! 


whore  -t  '  •  1  . 


A  two  point  (Jauss  integration  produces  from  (3.1)  the  approx¬ 
imate  - 


it 


e 


I 


1,-3  ..2  1  -1 ,  -1 . 2  1 ,  -12 

aril  r.w.  +,_r.  h  W.  +  rhr.  u. 

3  j  j  12  j  j  4  j  j 


-3  2^ 

+h  (hu.  +  wT)  r. 

J  J  J 


-hr  f ,w. 
2  j  J  J 


(3.3) 


for  which  the  values  of  u.,w.,u  ,w.  anil  w.  ari'  computed  rrom 
l  I  I  J  ;  .  I  J  I. 

u.  =  u,.'!  .  ,  w  =  u  ...  u .  =  ii  .  ,  w.  -  u  ■  .  ,  and  w.  *  u  ...  with 
J  'J  J  “  v  J  J  J  I  e  .1  J  f  j 

the  aid  of  the  numerical  vectors 


3l  -  ^(iWl.O.O/.jWl.O.O) 

J,  =  j‘(+2</3- 3,0, 0  t-'i/l , +2 v 3+ 1,0,0) 

-  1  .  -  <> 


-I  •,  =  r„(0, 9+4^3, 3+v'3,0,0,9+A^3,-3W3) 

I  1*-  lo 

’.‘J  2  -■  i(0,dV3,-uV3,O,O,3>/l.ll./)) 

Krom  (3.1)  we  produce 


(3.4) 


152 

l.  FR1KD 

II 

0) 

'V  „ 

;)ue 

i)2 irt, 
32uc 

-  =  'y  a.ii  !j)l+b  ■ 
J~j-J  J' 
j  =  l 

+d  .  ij> .  <p  j  +  i 

J-.l-j 

witli 

a .  = 
J 

2.-3 

3h  rj 

*  bJ 

=  rh  *  r  .  '  +  1 2h 

f>  .1 

-3  -2  J 

rjwj 

c.  =  4b  2r.w.  <1.  =2h  'r.  ,  e.  =  ^hr.  * 

J  J  J  J  J  I  2  j 


-2 


(3.7) 


J  J 


(3.8) 


To  assess  the  performance  of  our  element  we  use  it  to  com¬ 
pute  the  deflection  of  the  uniformly  loaded  (i.e.  f “const.), 
clamped  (i.e.  u(l)=w' (l)=0)  plate  with  an  immovable  (i.e.  u(l) 

=0)  edge.  For  f  *  10  we  compute,  with  Ne=2 , 3 , . . . , 7 ,  a  central 
deflection  w(0)  =  1 .  138993, 1  . 140754 ,  1.141714,  1.142070,  1.142220, 
1.142293;  meaning  a  relative  error  in  w(0)  equal  to  0. 068Ne~3 • 3 . 

A  Newton-Raphson  solution  of  the  nonlinear  stiffness  equation 
for  f  =  10  and  Nc  =  7  successively  comes  up  with  w(0)= 1 .874962 , 
1.394681  ,  1.182066,  1.143467,  1.142294,  1.142293,  liaving  star¬ 
ted  with  a  zero  deflection. 

2 

The  critical  thrust  for  the  plate  Is  given  by  pcr-  '  / 1 2  where 
3.39  is  the  first  root  of  Bessel's  function  equation  Ji('t)“0. 
Figure  7  traces  the  computed  central  deflection  of  a  simply  sup¬ 
ported  plate  bent  under  the  combined  action  13,101  of  a  lateral 
load  f  and  an  edge  thrust  p  that  exceeds  p,.r.  When  I  -  0  the 
unstable  trivial  solution  w  =  0  for  p  prr  is  absent  and  a  zero 
initial  deflection  can  be  chosen  for  the  Newton-Raphson  method. 
When  f  -  0  the  Newton-Raphson  method  must  start  with  a  nonzero 
initial  guess  but  proceeds  otherwise  as  bolorc  to  produce  the 
typical  bifurcation  curve  in  Fig.  7. 

Close  to  a  critical  point  at  which  K  '  is  nun  computable,  the 
condition  of  K  declines  and  the  Newton-Raphson  method  slows  down. 
It  is  our  experience,  though,  that  by  using  higher  precision 
computations  and  more  corrective  iterations  one  can  get  as  close 
as  it  is  only  numerically  meaningful  to  such  a  point. 


4.  RI1BBKR  MKMBRANK 

bet  the  generating  curve  of  the  undeformed  axisymmetric  mem¬ 
brane  be  given  In  the  (r,z)  plane  through  r=r(s)  and  z*z(s),  s 
being  the  arc  length.  Under  the  action  of  applied  forces  and 
prescribed  displacements  the  point  (r,z)  moves  to  the  deformed 
location  (x,y).  An  arc  element  ds  is  stretched  thereby  to  tin 
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FIG.  7.  Face  and  edge  forced  circular  plate. 

and  the  thickness  t  of  the  membrane  is  reduced  to  t. 


The  energy  density  function  of  the  membrane  is  ultimately  ex¬ 
pressed  in  terms  of  the  three  principal  stretch  ratios 


.  -  «L»  .  \  ,  ?*x  ,  =  t 


1  '  ds"  ’  2  ‘  Hr  ’  3  ”  t 

and  if  the  deformation  is  incompressible 
ds  =  (dx2+dy2)'i  the  stretch  ratios  become 


(4.1) 


With 


(x 


>2  ♦  v»-’)2  v  -  x  i  ,  .  . 

12 


(4.2) 


where  (  )'  -  d/ds. 


Asstiming  a  membrane  made  of  Mooney  material  its  total  poten¬ 
tial  energy  acquires  the  lorni 

tr(x,y)  =  2irut  J  I  ( 1  -3)  t  >(  I  ,(-3)  !rds+i  ^  I  x“y'ds  (4.3) 
0  1  "  -  ,;l  t) 


where  p  is  the  pressure,  :i  and  <  material  constants,  and  where 
1|  and  I,(  are  the  strain  invariants 


2.2  2 

1 1  =  ^  1  +  2  +  ‘l 

2  2  2  2  2 

*  V.  ' 


(4.4) 


with  >  .j  =  From  here  on  we  replace  p/|it  by  |*  am!  assume 

tint l  2* in.  *  I. 


l-AKCi:  I.I.AS  I  l<  Dl.idKMM  lo-.:; 


b  = 


~rds  h 


ri)s  1, 


>  r.^'.vV 

i  I  ’  ''  '' 


/  , 

i  ;l 


J  -1! 

r .  x . 

I  ! 


X  v'.l: 


A  Vi 

i  I 


I'fpr.it  ii!  cl  1 1  !  t  T«nl  i .  i  f  ion  o!  r%.  with  rc-spri-t  t«»  u  tui’iiislu 


7‘  u. 

•u,. 


z 


,]=i 


a  :  H  I)  4 1  • 

II  i  i  i  i 


with 


_'I  •> 


=  2v  i°-  n  Vn+  v 


=  -r.y  .  (1  Cl  *  )  4 

I  I  .1  »  I  2  I  2  I  j 


C  .  = 

-J 


2,,'2JM-'U  +  h',xA. 


r  i 


;ind 


■i  :  i  !.  i  . 

i  i  i  i  i  i  iii 

hi 


)  )-J  ,i-.i  .]  -j  j  -  -j' 


h  f  .! .  ) 
i  i  J  i  i 


■< .  ■>  > 


(4. 10) 


(4.11) 


with 
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i.  r:ii:i> 


a .  -  i'ti  '  r  .  (1  - 


••  '  .  'tv 

j---  J  ) 
2  6 

•’  i  i  i 

x :  ■  ’  - 


h.  21.  *  r .  ( t  -  -  -  -  --> 

I  I  ? 


)  ( I  ♦  •  .  I 


>n  + 

2  i 


•  -  2hr  1  (1+  1  \  “  '*)  <  1  +  •;  .  )  4  !ipv '  . 

I  i  I  i  i  I  i  i 


a.  4x.  •,.(  •  + 

]  J  ij  I  |  2) 


f4.1  2) 


i  '  i  ’  i  I  i  i  "i 


i .  =  hi,  n  .  •  • .  > 

i  i  i  i  1 1  -  i  .  i 


()iir  first  appi irat  inn  of  the  rubber  membrane  clement  is  to 
compute  the  inflated  shape  1  1  i  of  a  unit  disc  (  •  =  0.1)  lor  a 
pressure  p  *  5  and  an  inplane  stretching  of  tin  edge  to 
xfl)  ■  1.1.  The  disc  is  suh  st  ant  i  a  1  I  v  dclormcd  hy  t  hi?;  high 
pressure,  and  for  a  uniform  layout  of  N'c  I  mile  elements  we  com¬ 
pute,  correspond  iup.  to  fie  -  1,2,1,4,11),  a  polar  rise  y  (0)*  1  .2546, 
I  .4278,  1.42')'),  1. 4104,  1.4104.  We  reach  this  last  value  of 
y(0)  with  the1  Newton-Kaphson  scheme  tliat  successively  computes 
y (0)=2 . 31 875 ,  1.44782,  1  .44302,  1  .4  10V),  1.4  1040,  1.4  3040. 

Figure  8  shows  inflated  .hopes  ot  the  di  c  lor  a  pressure 
p  =  1,2,.  ...7.  Correspond  ing  to  these  pressures  arc-  the  polar 
stretch  ratios  >„  =  1.144,  1.233,  1.171,  1.648,  2.430,  4.825, 
6.552.  The  global  stiffness  matrix  K  is  found  to  he  positive 
definite  for  the  shapes  in  Fin.  8  and  we  conclude  that  the  disc 
is  in  stable  equilibrium. 

Figure  9  shows  a  torus  1 6  I  (  •  ~  0.25),  generated  bv  r  “  2  + 
cos  s,  7.  -  sin  s  0  '  s  " ,  inflated  by  a  pressure  that  in¬ 
creases  in  ten  equal  steps  from  p  -  0  to  a  critical  p  =  2.185. 

As  the  pressure  approaches  this  last  value  of  p  the  lowest  ei¬ 
genvalue  of  the  global  stiffness  matrix  K  nears  zero  Indicating 
a  decline  in  stability. 


i  jk  i.  ikikp 

'■'inure  IP  shows  hulninn  '  f>J  o;  a  tutu-  (  =  0)  ,  stretched 

both  axially  and  c  ircurr.f  crmt  1  y ,  and  inflated.  The  curves 
in  Fig.  10  are  for  a  pressure  that  increases  in  ten  equal  steps 
from  p  =  0  to  the  erit  icnl  p  -  O.UT. 


/"'///W 

/  /  '/  ////// 

'  //sy/M 


!••[(;.  1!>.  Itul.’.im;  of  stretched  and  in'lated  tube 
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Discrete  Gauss  integration  of  the  element  total  potential  energy  is  applied  to  the  formation  of  a 
cubic-cubic  C'  extensible  completely  nonlinear  curved  beam  finite  element.  The  versatility,  accuracy, 
effectiveness,  and  robustness  of  the  element,  and  the  Newton-Raphson  technique  used  to  solve  the 
nonlinear  algebraic  stiffness  equation  set  up  with  it  is  numerically  demonstrated  by  computations  of  the 
nonlinear  equilibrium  stability  and  motion  of  beams  and  rings. 


1.  Introduction 

A  nonlinear  finite  element  discretization  technique  based  on  the  approximate  Gauss 
integration  of  a  nonquadratic  energy  density  function,  successful  in  the  nonlinear  computation 
[1]  of  the  straight  and  curved  inextensible  elastica,  the  circular  plate,  and  the  axisymmetric 
rubber  membrane,  is  applied  here  to  nonlinear  equilibrium  stability  and  motion  analysis  of  the 
extensible  curved  beam. 

A  cubic-cubic  C*  element  is  developed  in  detail  (a  numerically  integrated  cubic-cubic  beam 
element  for  large  displacements  is  available  in  the  MARC  program  [2])  and  is  computationally 
tested  for  accuracy  and  effectiveness  on  the  particular  large  displacement  problems  of  a  tip 
loaded  straight  beam,  a  closed  ring  compressed  by  two  equal  and  opposite  forces,  a  circular 
ring  under  post  critical  hydrostatic  pressure,  and  the  large  amplitude  vibrations  of  free  and 
fixed  beams  and  rings. 


2.  Finite  dement 


With  reference  to  Fig.  1  let  a  point  on  the  deflected  beam  by  marked  by  (x,  y),  x  =  x(s)  and 
y  -  y(s),  s  being  the  distance  measure  along  the  original  curved  beam.  Let  further  e  denote 
the  axial  strain  and  k  the  curvature  of  deflected  beam.  Then  in  terms  of  x  and  y 


and 


e  =  (x,2+  y’2)m-  1 

Y*  V*  —  V*JC* 


(1) 

(2) 
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Fig.  1.  Geometry  and  element  of  a  largely  bent  beam. 

with  the  usual  notation  x'  =  dx/d s  and  x"  -  d2x/ds2.  The  initial  curvature  of  the  beam  is 
expressed  in  terms  of  the  original  stope  6  as  k0  =  dfl/d s. 

A  bent  and  stretched  elastic  beam  of  length  /  that  is  under  the  action  of  distributed  forces  / 
and  g  in  the  x  and  y  direction,  respectively,  is  in  possession  of  a  total  potential  energy 

-rr(x,  y)  =  j  El  f  (k-  Kofds  +  \EA  f  e2ds-  f  (fx  +  gy)ds  (3) 

Jo  Jo  Jo 

or  when  the  beam  is  of  thickness  t  and  unit  width 

tt{x,  y)  =  nEt3 [l  Jg  (*  ~  Kofds  +  |c  £  e2  ds J  -  (fx  +  gy)ds  (4) 

where  c  =  12/ 12.  Hence  forward  we  shall  assume,  for  typographical  briefness,  that  Et3 f  12  =  1. 

For  the  approximation  of  ir(x,  y)  in  (4)  we  propose  a  C1  cubic-cubic  finite  element  and  a 
three  Gauss  point  integration,  which  appears  to  be  the  minimal  integration  scheme  to  maintain 
the  full  element  accuracy  inherent  in  the  cubic  interpolation,  while  averting  spurious  zero 
energy  modes.  We  shall  also  consider  in  the  paper  the  possibility  of  integrating  the  axial  strain 
energy  part  of  ir(x,  y)  with  only  two  Gauss  points,  but  it  appears  that  an  all  out  three  point 
integration  of  ir(x,  y)  is  preferred. 

To  prepare  for  the  numerical  integration  the  typical  finite  element  is  mapped  from  s  to  f  by 
s  =  Si  +  hi;,  0  s  £  s  1,  so  that  ds  =  h  d|  and  (  )'  =  h~’(  )  and  (  y  =  h~'{ "),  where  the  dot  means 
d/d£  From  the  element  nodal  values  vector  (see  Fig.  1) 


u'a  =  (xu  x,,  yu  yu  x2,  x2,  y2,  yi) ;  (5) 

x  and  y  are  both  cubically  interpolated  inside  each  element  by 

x  *  and  y  -  u\ift,  (6) 

<t>  and  0  being  the  shape  function  vectors 

<t>'  =  (4>u  <t>2, 0, 0,  fa  <£4, 0, 0)  and  ift'  =  (0, 0,  <f»u  fa,  0, 0,  <fc,  <^4)  (7) 
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in  which 


<a,=i-3*2+2£3,  <fe=f-2f+i3,  <*>3=3f2-2i3,  <t>A=-e+e. 


Since  the  typical  element  extends  between  £  =  0  and  (=1,  and  since  ds  =  h  d£  the  ap¬ 
proximate  integration  of  the  element  total  potential  energy  ire(x,  y )  is  of  the  form 

TTc  =  h  2  wAHkj  -  Kojf  +  \ce)  -  fix,  -  gi.v,]  (9) 

where  here,  specifically, 

*Vi  -  Wj  =  "fl,  W2  =  ^  (10) 

with  the  three  Gauss  points  GIt  G2,  G3  being  at 

f,.3  —  w(5  +  VT5)  and  £2  =  \  (11) 

the  upper  sign  of  vT5  referring  to  j  =  1  and  the  lower  sign  of  V15  referring  to  j  =  3. 

The  integration  point  values  of  the  curvature  k,,  the  strain  e,,  and  the  coordinates  x,  and  y, 
needed  in  irc  in  (9)  are  obtained  from  the  nodal  values  vector  ue  through 

X,  =  x,  =  x,  =  u% 

y(  =  y,  =  «'.*>,  if)  =  Uttii 

where  <f»,  stands  briefly  for  *(£),  etc.  Here,  for  the  three-point  Gauss  integration  we  have 
from  (6),  (7),  (8)  and  (11)  that 

4>\j  «  i5o(50  ±  12VT5,  5  ±  Vl5,  0,  0,  50  +  I2VT5,  -5  ±  Vl5,  0,  0) , 

*2  =  1(4,  1,0,0,  4,  -1,0,  0), 

* u  =  ^(-6,  2 ±  Vl5,  0,  0,  6,  2  *  Vl5,  0,  0), 

**2  =  i(-6,  - 1, 0, 0, 6,  - 1, 0, 0) ,  (i  3  ‘ 

*',.3  =  l(*6Vl5,  -5  ±  3Vl5,  0,  0,  ±6VT5,  5  *  3Vl5,  0,  0) , 

*‘2  =  (0,-1,  0,0,0,  1,0,0) 

and  ___ 

*u  =  i»(0,  0,  50  ±  12V15,  5  ±  Vl5,  0,  0,  505: 12VI5,  -5  ±  VB) , 

*’z  =  i(0,  0,4,  1,0,0,  4,-1), 

*u  =  w(0,  0,  -6,  2  ±  VT5,  0,  0,  6,  2  ¥  Vl5) ,  (14) 

*i»i(0,  0,-6, -1,0,0,  6,-1), 

*’,.3  «  j(0,  0,  *  6Vl5,  -5 ± 3VT5,  0,  0,  ±6Vl5,  5  +  3Vl5) , 

*2  =  (0,  0,  0,  -1,0,  0,0,1). 
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Differentiation  of  ire  in  (9)  with  respect  to  the  element  vector  of  nodal  values  ue  creates  the 
element  gradient  vector  ge.  Further  differentiation  of  ge  with  respect  to  ue  generates  the 
nonlinear  element  stiffness  matrix  kc.  Before  doing  that  we  shall  introduce  some  notational 
abridgements.  First  we  rewrite  the  axial  strain  €  and  the  curvature  k  of  the  beam  as 

e  =  1  and  k  =  a/3-3/2  (15) 


with  (in  fact  for  a  nearly  inextensible  beam  we  may  set  /3  =  1  in  k) 
a  =  xy-  xy  and  /3  =  x2  +  y2  . 


(16) 


Secondly,  we  wish  to  employ  in  this  section  a  prime  to  denote  differentiation  with  respect  to  the 
vector  uc  so  that  e’  and  k’  are  vectors,  and  e",  k",  eV  and  kV'  are  matrices.  Now,  since 
x)  =  4>,  and  y)  =  ^ 

3 

<  =  ge  =  h  X  W/U*,  “  «o/)k'  +  ceje'j  -  fa  -  gjift,] 


/-> 


in  which 


and 


a',  =  x,$,  +  thy,  -  x,4fi  -  4>,y, ,  ft  =  2(x^  +  y,4>t) > 

*'i= 


= \h'iK'npi . 

Next  we  find  that 

7r"  =  lct  =  ll  2  -  *0))*”  +  K’,K?  +  C^"  +  €•€■*) 

/-I 

where,  in  terms  of  and  ft, 

and 

k?  =  Pi7a[fia’;-32«AP}-to(fiW  +  aP7)  +  ¥«*ft0j] , 

both  computed  from 

«7  =  •M')  +  ^5  -  M)  ~ 


and 


(17) 

(18) 

(19) 

(20) 

(21) 

(22) 

(23) 

(24) 

(25) 


which  are  numerical  symmetric  matrices  of  dimension  8x8. 

Assembly  of  g ,  and  km  into  the  global  g  and  K  is  the  same  here  as  for  linear  finite  elements 
except  that  an  iterative  method,  say  Newton-Raphson,  is  needed  to  solve  the  nonlinear 
discrete  equation  of  equilibrium  g  *  0.  With  the  Newton-Raphson  method  an  initial  guess  Uo 
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is  improved  into 

u,  =  Uo-  Ko'go,  (26) 

etc.,  until  convergence. 


3.  Tests  for  effectiveness 

First  to  concern  us  is  the  accuracy  of  the  element,  occurrence  of  spurious  modes  devoid  of 
energy,  and  the  performance  of  the  Newton-Raphson  method  in  the  presence  of  a  large  axial 
elastic  constant  c.  We  perform  tests  in  this  respect  with  two  numerical  integration  schemes: 
one  that  integrates  the  bending  part  of  the  element  total  potential  energy  with  three  Gauss 
points,  and  the  stretching  part  of  irt  with  two  points  (this  will  be  referred  to  as  the  3-2  scheme); 
and  another  procedure  that  integrates  the  entire  irc  with  three  Gauss  points  (a  3-3  scheme). 
To  start,  we  compute  the  deflection  of  an  originally  straight  unit  cantilever  beam  with  a  tip 
force  P  -5,  shown  in  Fig.  3,  and  imposed  boundary  conditions  jt(0)  =  y(0)  =  y'(0)  =  0. 
Standard  assembly  of  all  the  element  stiffness  matrices  produces  a  global,  displacement 
dependent,  stiffness  matrix  that  can  be  put  in  the  form 

K=  Kb  +  cK%  (27) 

‘  with  Kb  and  K,  constituting  the  bending  and  stretching  parts  of  K,  respectively. 

All  subsequent  computations  are  carried  out  with  7  and  14  element  discretizations  that  give 
rise  to  27  and  57  degrees  of  freedom,  correspondingly.  The  3-2  integration  scheme  is  found  to 
•  produce  a  K,  matrix  with  15  zero  eigenvalues  for  a  7  element  discretization,  and  29  zero 
eigenvalues  for  the  14  element  discretization.  Imposition  of  the  boundary  conditions  x(0)  = 
y(0)  =  y'(0)  still  leaves  the  total  global  stiffness  matrix  K  with  one  spurious  zero  eigenvalue  in 
its  straight  configuration,  which  disappears  with  bending,  and  which  does  not  seem  to  heap  any 
difficulties  upon  the  working  of  the  Newton-Raphson  method  even  with  a  straight  (x  =  s,y  = 
0)  initial  guess.  Table  1  lists  the  largest  eigenvalue  A*  of  the  global  stiffness  matrix  K  as  it 
varies  with  the  number  of  elements  Ne  and  the  axial  elastic  constant  c.  Table  2  lists  the  tip 
coordinates  x(l)  and  y(l)  computed  (with  some  16  significant  digits)  with  the  Newton- 
Raphson  method  for  different  values  of  c. 

Raising  the  integration  scheme  to  a  3-3  level  produces  an  element  stiffness  matrix  kt  that 
assembles  into  a  global  K  with  a  K,  part  that  has  14  zero  eigenvalues  for  N,  -  7  and  28  zero 
eigenvalues  for  N.  =  14.  No  spurious  zero  eigenvalues  occur  anymore  in  the  assembled  and 
constrained  K ,  the  extremal  eigenvalues  of  which  are  listed  in  Table  3  in  their  dependence 
upon  N,  and  c.  Table  4  lists  the  computed  tip  coordinates  x(l)  and  y(l)  of  the  originally 
straight  beam  discretized  with  3-3  integration  elements  and  tip  loaded  with  /*  =  5,  as  they 
become  improved  with  the  Newton-Raphson  method,  for  different  values  of  the  axial  elastic 
constant  c,  and  a  different  number  of  finite  elements  Nm. 

A  remarkably  propitious  conclusion  emerges  from  Tables  2  and  4:  that  the  Newton- 
Raphson  method  is  only  slightly  affected  by  the  large  values  of  c. 
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Table  1 

Largest  eigenvalue  A  S  of  K  formed  with  3-2  integration  elements 
in  its  dependence  upon  the  number  of  elements  N.  and  the  axial 
elastic  constant  c 


I'¬ 

ll 

2? 

ll 

£ 

c=  103 

c  =  104 

c  =  105 

c  =  105 

AS 

0.27  105 

0.27  106 

0.27  10’ 

0.55  !07 

Table  2 

Convergence  of  the  tip  coordinates  x(l)  and  y(l)  with  the  Newton-Raphson  (NR)  iterative  cycles.  Straight  beam  tip 
loaded  with  P  =  5  and  discretized  with  3-2  integration  elements 


N. 

,  =  7 

N„ 

=  14 

c 

=  103 

c  = 

■  104 

c 

=  105 

c  = 

=  105 

x(l) 

yd) 

*(1) 

y(D 

x(\) 

yd) 

x(l) 

yd) 

1 

1.0000000 

1.6830084 

1.0000000 

1.6830084 

1.0000000 

1.6830084 

1.0000000 

1.6924615 

2 

0.55703212 

0.77655095 

0.55205671 

0.77962496 

0.55134776 

0.78051978 

0.54965834 

0.78171462 

3 

0.61384381 

0.72384822 

0.61361089 

0.72184756 

0.61276418 

0.72211633 

0.61403859 

0.72174536 

4 

0.61061737 

0.71912615 

0.60820549 

0.71735096 

0.60713127 

0.71758132 

0.60799521 

0.71731317 

5 

0.61180749 

0.71817148 

0.61099979 

0.71498461 

0.61090756 

0.71467798 

0.61042256 

0.71501707 

6 

0.61180735 

0.71816995 

0.61099710 

0.71497673 

0.61091337 

0.71465933 

0.61041886 

0.71501085 

7 

0.61180735 

0.71816995 

0.61099714 

0.71497671 

0.61091606 

0.71465766 

0.61041902 

0.71501073 

8 

0.61099714 

0.71497671 

0.61091606 

0.71465766 

0.61041902 

0.71501073 

Table  3 

Lowest  (1st)  and  highest  (N  th)  eigenvalues  A  *  and  A  S  of  the  assembled  global  stiffness  matrix 
K  for  the  straight  beam  discretized  with  N.  3-3  integration  elements 


N.  =  7 

N.«  14 

c-  10* 

c=  10* 

c-  105 

c-  10s 

Af  Aft 

Af  Aft 

Af  AS 

Af 

AS 

1.25  0.32 10* 

1.25  0.32 10“ 

1.25  0.32 107 

0.75 

0.66 107 

/ 
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Table  4 

Convergence  of  the  tip  coordinates  x(l)  and  y(l)  with  the  Newton-Raphson  (NR)  iterative  cycles.  Straight  beam  tip 
loaded  with  P  =  5,  and  discretized  with  3-3  integration  elements 


N. 

=  7 

N. 

=  14 

c 

=  103 

c  = 

=  10* 

c 

=  10s 

c  = 

=  105 

•*0) 

y(l) 

*(1) 

yd) 

x(l) 

yd) 

*0) 

yd) 

1 

1.0000000 

1.6651444 

1.0000000 

1.6651444 

1.0000000 

1.6651444 

1. 0000000 

1.6636502 

2 

0.56022070 

0.77453293 

0.55526060 

0.77770589 

0.55452073 

0.77870079 

0.55478261 

0.77847397 

3 

0.61529505 

0.72278893 

0.61557099 

0.72053522 

0.61541080 

0.72053192 

0.61594304 

0.72015670 

4 

0.61221822 

0.71820317 

0.61007604 

0.71628612 

0.60969777 

0.71624185 

0.61019765 

0.71596357 

5 

0.61335240 

0.71729088 

0.61260781 

0.71402538 

0.61297276 

0.71312999 

0.61251874 

0.71374747 

6 

0.61335228 

0.71728948 

0.61260543 

0.71401851 

0.61296877 

0.71311681 

0.61261565 

0.71374152 

7 

0.61335228 

0.71728948 

0.61260545 

0.71401850 

0.61296970 

0.71311611 

0.61251582 

0.71374140 

8 

0.61260545 

0.71401850 

0.61296970 

0.71311611 

0.61251582 

0.71374140 

9 

0.61237501 

0.71378936 

The  (stable)  equilibrium  configuration  computed  for  the  tip  loaded  beam  with  the  3-2 
integration  scheme  is  shown  in  Fig.  2.  Fig.  3  shows  the  computed  curvature  distribution  k(s) 
for  the  same  beam,  and  Fig.  4  shows  the  tension  p  along  this  beam  computed  from  equilibrium 
considerations  (p  =  Py')  and  from  Hook’s  law  (p  =  ce).  We  remark  in  Fig.  3  that  even  though 
the  elastically  computed  tension  oscillates  violently  inside  each  element,  at  the  two  Gauss 
points  Gi  and  G2,  at  which  the  stretching  energy  is  sampled,  the  statically  and  elastically 


computed  tensions  agree  and  are  accurate  to 
requirement. 


Fig.  2.  Deep  bending  of  a  tip  loaded  beam  computed 
with  3-2  integration  elements. 


a  degree  that  is  sure  to  satisfy  any  practical 


Fig.  3.  Curvature  distribution  for  the  beam  of  Fig.  2. 


Fig.  7.  As  in  Fig.  6  but  with  c  =  10*. 


3-3  Integration 


k-( x'y“-y'x')/(x‘x,*yy') 


c-IOOO 


Fig.  8.  Curvature  distribution  for  the  beam  of  Fig.  2  discretized  with  3-3  integration  elements. 

Figs.  5, 6  and  7  trace  the  computed  tensions  for  the  tip  loaded  beam  discretized  with  3-3 
integration  elements.  Fig.  8  shows  the,  here  smoother  looking,  computed  curvature  dis¬ 
tribution  for  the  same  beam. 


4.  Pressed  ring 

The  elastic  deflection  of  a  thin  inextensible  circular  ring  pressed  by  two  equal  and  opposite 
forces  has  been  previously  computed  [1}  with  quadratic  C°  finite  elements.  In  this  section  we 
present  similar  results  obtained  with  the  cubic-cubic  C 1  element.  Fig.  9  shows  the  equilibrium 
configurations  that  the  ring  rssum^s  when  pressed  by  a  force  P  that  increases  in  steps  of  2 
from  0  to  38.  Fig.  10  follows  the  '  Ming  of  the  gap  along  the  diameter  of  the  ring  with  the 
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Fig.  10.  Computed  sinking  of  ring's  top  point  with  the 
Fig.  9.  Point  pressed  circular  ring.  top  force  P. 

increase  of  P.  We  see  in  Fig.  10  that  our  extensible  element  has  a  slight  extra  flexibility  to  it; 
contact  of  the  pressed  points  is  reached  with  a  force  P  that  is  somewhat  less  than  the 
theoretical  [3-5]  P  =  4-rr2  for  a  ring  of  radius  1/ir. 

5.  Ring  under  hydrostatic  pressure 

To  account  for  the  action  of  an  external  uniform  pressure  (nonuniform  is  actually  as  easy)  q 
we  have  to  add  to  rr(x,  y)  in  (3)  the  pressure's  potential 

7r(x,  y)  =  q  f  xy'  ds  (28) 

JO 

or  approximately,  for  the  typical  element, 

w wy,  (29) 

that  leads  to 

g.  =  it;  =  q  X  W/Ww  +  (30) 

and  then  to 

fc«  =  irnt  =  q  £  +  fo))  (31) 

i- 1 

which  are  added  to  gt  and  kt  in  (17)  and  (21). 
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Fig.  11.  Deformation  of  a  circular  ring  under  the  action  Fig.  12.  Displacement  of  the  ring's  top  point  under  the 

of  an  external  pressure  q.  action  of  an  external  pressure  q  and  a  top  force  P. 

The  critical  pressure  on  a  dosed  circular  inextensible  ring  is  at  q„  =  3/r2,  or  in  our  case 
where  r  =  1/tt,  qct  =  3ir2  =  29.609.  Fig.  1 1  shows  the  collapse  of  the  rir.g,  discretized  with  seven 
elements,  under  the  action  of  an  increasing,  as  listed,  pressure.  Fig.  12  compares  the  variation 
of  the  maximal  deflection  w  caused  by  the  external  hydrostatic  pressure  q  and  a  point  force  P, 
computed  with  our  element,  with  results  obtained  from  an  analytic  quadratic  theory  [6-9].  A 
high  c  value  is  chosen  here  to  make  the  comparison  with  the  inextensible  theory  more 
meaningful. 


6.  Vibrating  beams 

If  g(u)  denotes,  as  before,  the  global  gradient  of  the  total  potential  energy  and  M  the 
beam’s  global  mass  matrix,  then  its  equation  of  motion  is  written  as 

g(u)  +  Mu  =  0  (32) 

where  u  denotes  velocity  and  ii  acceleration,  that  we  propose  to  numerically  solve  with  the 
Newmark  scheme  [10] 

u,  -  Uo+  TUo  +  {t2Uo,  u,  =  Mo  +  5T(t«o+  Ml)  (33) 

for  a  time  step  size  t.  By  virtue  of  (32),  (33)  becomes 

M,  =  Mo+  TWo-5T2Af"'go,  Ml  =  «o~  iTA/-'(go+  g,).  (34) 

In  our  subsequent  dynamic  computations  we  shall  exclusively  use  the  consistent  [10]  beam 
mass  matrix. 
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We  can  foresee  the  possibility  of  a  simply  supported  beam  with  free  axial  motion  executing 
very  large  amplitude  vibrations  even  when  nearly  inextensional.  To  numerically  study  such 
motions  we  impart  the  originally  straight  (x  =  s,  y  =  0)  beam  an  initial  velocity 

xo  =  0,  yo  =  ~  sin  -its  (35) 

with  magnitude  determined  by  the  factor  a.  Fig.  13  shows  the  appearance  of  the  beam  at  time 
intervals  of  m  as  it  hurls  up  towards  its  ultimate  flexed  position.  In  Fig.  14  the  beam's  central 
elevation  is  traced  against  time  for  different  values  of  a  in  (35).  When  the  vibration  amplitude 
is  small,  say  for  a  <  2,  the  computed  time  it  takes  the  crest  of  the  beam  to  rise  and  fall  back  to 
zero  is  close  to  the  theoretical  half  period  value  of  I/tt.  A  growing  amplitude  is  predicted  by 
our  computations  to  cause  a  hardening  of  the  beam  and  a  shorter  period  of  vibration.  This 
does  not  sit  well  with  other  published  results  [11, 12],  but  meaningful  comparisons  are  anyway 
not  easy  here.  First,  the  approximations  of  the  analytic  approaches  are  sensitive  [11]  to  very 
large  displacements;  and  also,  one  must  bear  in  mind  that  for  a  periodic  solution  to  exist  in  the 
nonlinear  range,  particular,  not  easy  to  come  by,  initial  conditions  must  be  at  hand. 

To  study  the  large  movements  of  a  free-free  beam  we  set  it  in  motion  with  the  initial 
velocities 


jco  =  0  and  y0  =  a[cos  As  +  cosh  As  +  a(sin  As  +  sinh  As) 


where  A  =  4.7300408  and  where 


sin  |A  -  sin 
cos  j  A  +  cosh  i  A 


-1.018. 


(36) 

(37) 


Fig.  13.  Movement  of  a  simply  supported  axially  un¬ 
constrained  beam. 


Fig.  14.  Movement  of  the  beam's  central  point  for 
different  magnitudes  of  the  initial  velocity  in  equation 

(35). 
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Fig.  15.  Movement  of  a  free-free  beam.  Fig.  16.  Progress  of  the  total  bend  in  a  free-free  beam 

for  different  values  of  a  in  equation  (36). 


Fig.  15  shows  a  sequence  of  still  pictures  computationally  taken  of  the  beam  between  its  first 
and  second  rest  positions  for  a  velocity  factor  a  =  2.  Fig.  16  traces  the  total  bend  y(0)-  y(j)  of 
the  beam  as  it  progresses  with  time.  Again,  for  small  amplitudes  the  free-free  beam’s 
computed  half  period  is  near  ir/A2  predicted  by  the  inextensible  small  displacements  theory, 
but  as  a,  and  with  it  the  amplitude,  keeps  growing  the  beam  becomes  softer  with  longer 
periods,  in  agreement  with  the  computations  of  Takahashi  [13]  but  opposite  those  of  Wagner 
[14]. 


7.  Vibrating  ring 

For  the  circular  ring  it  is  more  convenient  to  deal  with  the  normal  displacement  w  and  the 
tangential  displacement  v,  that  are  related  to  x  and  y  (see  Fig.  17)  by 


av 
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Fig.  17.  Vibrating  ring.  Fig.  18.  Movement  of  w(0)  of  the  ring  in  Fig.  17. 

or  in  terms  of  the  Cartesian  coordinates 

io  =  tvocos  0  -  t)osin  6,  y»  =  WoSin  6  +  0()cos  0  (41 ) 

with  0  s  9  <  jit. 

Fig.  17  shows  a  train  of  snapshots  computationally  taken  of  the  beam  at  m  time  intervals. 
Fig.  18  traces  w(0)  with  time  for  different  magnitudes  of  the  initial  velocity.  When  a,  and  with 
it  the  displacement  amplitude,  is  small  the  ring  executes  half  a  period  in  a  time  nearly  equal  to 
t^e  theoretical  inextensible  V5ftm,  but  as  the  amplitude  grows  a  slight  softening  is  com¬ 
putationally  detected  for  the  ring  in  agreement  with  Evensen’s  [15]  reporting.  Fig.  19  shows 
the  influence  of  c  on  the  motion  of  the  ring  started  with  h>„  and  v0  in  (40). 


Fig.  19.  Influence  of  c  on  the  movement  of  w(0)  for  the  ring  in  Fig.  17. 
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8.  Nonconservative  loads 

Beck’s  problem  of  a  slender  cantilever  loaded  by  tip  follower  forces  [16]  is  perhaps  of  little 
practical  value  as  its  real  occurrence  is  somewhat  farfetched  but  it  will  serve  as  a  good 
example  to  describe  the  application  of  finite  elements  to  problems  with  no  potential  for  the 
load.  Fig.  20  shows  the  arrangement  of  the  tip  forces  that  consist  of  a  tangential  follower  force 
P  and  a  normal  follower  force  Q.  When  the  load  cannot  be  derived  from  a  potential,  as  here, 
we  are  forced  to  forego  the  addition  of  appropriate  work  terms  to  the  total  potential  energy 
and  must  introduce  instead  the  forces  directly  into  the  gradient  and  consequently  into  the 
(nonlinear)  equations  of  equilibrium.  For  our  present  beam  element  this  operation  is 
straightforward.  Indeed,  since  the  gradient  entries  that  correspond  to  the  x  and  y  nodal  values 
express  the  force  sums  in  these  direction,  all  we  have  to  do  here  is  to  add  the  negative 
horizontal 

H  =  P  cos  6  +  O  sin  0  (42) 

component,  and  the  negative  vertical 

V  =  P  sin  0  -  O  cos  0  (43) 

component  of  the  tip  forces  to  the  entries  of  g  that  correspond  to  the  tip  x  and  y,  say  gN-3  and 
gs  i  if  the  tip  node  is  the  last. 

When  we  know  the  beam  to  be  nearly  inextensible  the  substitutions 

sin  0  =  /T'y(l)  and  cos  0  =  /r'jr(l)  (44) 

can  be  made,  with  which  H  and  V  become 

H  =  h  '(PuN-2+  Qun)  and  V  =  h~'(Pus  -  Ous-2)  (45) 


fig.  20.  Cantilever  bent  by  large  follower  forces. 


fig.  21.  Same  as  fig.  20  but  with  P  m  21. 
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that  we  add  to  gN~3  and  gN-u  respectively.  To  prepare  for  the  Newton-Raphson  solution  of 
g(u)  -  0  the  four  terms 


dUs-2 


h'Q , 


dgs-i 

9Us-i 


h'Q, 


dgN-i 

duN 


h'P 


(46) 


have  to  be  added,  nonsymmetrically,  to  the  global  stiffness  matrix  K  at  the  addresses  (N  -  3, 
N  -  3),  (N  -  3,  N),  (N  -  1,  N  -  3)  and  (N  -  1,  N),  respectively.  Fig.  20  and  21  show  computed 
equilibrium  configurations  for  the  elastic  column  bent  by  the  follower  forces  P  and  Q.  No 
nontrivial  equilibrium  configuration  is  found  for  0  =  0. 
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samerical  iatograiioa  it  a  ride  issue  to  this  to 
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I  with  ank  wfiasmtat  ■  showi  to  be  aa  cf  active  ami 
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aad  that  (high  order)  exact  fcugratioa  of  togb  order  fiaitc 
viriatioa  principle  of  ctoetical  etosticity,  metoding  a  biph 


l  nrruooucnoN  greater  detail  the,  so  often  misquoted,  arguments  of 

to  1*74  the  author  wrote  a  paper  [1]  on  the  finite  element  Ref.[l]  and  to  dearly  point  out  that  unlike  the  penalty 

mriysis  of  nearly  incompressible  elastic  solids  in  which  methods  that  rely  on  huge,  machine  dependent,  bulk 

le  advocated  the  gradual  introduction  of  incompres-  modulus  and  reduced  integration  that  supposedly  release 

ribity  by  way  of  a  bulk  modulus  that  is  allowed  to  volume  preserving  modes,  the  present  technique  does 

■crease  only  when  accompanied  by  a  proper  mesh  not  need  them. 

■faction  in  such  a  way  as  to  balance  the  discretization  ____ 

■tors  with  the  compressibility  errors.  Progressive  «•  1  ■tCGumsHBUSHtiai 

rr*Tn  of  mcompressMity  is  called  for  in  the  finite  b  order  to  make  our  arguments  as  explicit  as  possMe 


'  gaoitm  of  incompressibility  is  called  for  in  be  finite  b  order  to  make  our  arguments  as  explicit  as  possMe 
dnaeat  displacement  method  because  for  any  given  we  shall  pay  close  attention  to  the  simple  problem  of  an 

■sth  a  too  large  bulk  modulus  can  either  cause  a  dis-  elastic  unit  hollow  sphere,  as  m  Fig.  1  stressed  by  a  unit 

onerously  iD  conditioned  stiffness  equation  or  an  overly  internal  pressure,  but  free  of  stress  on  its  outer  skin. 

Miff  (locked)  finite  element  model.  Exact  enforcement  of  Here,  when  the  shear  modulus  equals  one 

■compressibility  via  a  minimum  variational  principle  aad 

Hynomisl  shape  functions  is  often  possible  only  at  the  I  .  1“^  r1  .■  ■*  m 

hid  state  of  zero  displacements.  And  if  the  com*  '  *  '  *1.4?  2(1  +  •>)  r  l-o’  ' 

PMWiuual  model  is  tough  why  should  one  use  anyway  a 

k  hrge  bulk  modulus  and  disproportiooatly  smaO  coos-  o  being  the  internal  radius  and  v  the  Poisson  ratio.  With 
.  pess&Qity  errors? 

The  paper,  which  included  computations  of  aa  elastic  ....  1» 

.  tohere  came  upon  the  additiooal  observation  that  l-2r'  11 

■■erical  integratioa.  not  necessarily  reduced,  can 

-  "mm*  otter  locking.  - Equation  (1)  becomes _ 

;.  ,  atoee  then  the  subject  of  finite  dement  analysis  of 
■MtopretiiWe  materials  both  solid[2)  and  fcrid  [3>-7]  has  .  .  p 

kp<  central  stage  in  computational  continuum  mechanics  ■<r.z)-u(r.-)  +  5£5r.  (3) 


■compressible  materials  both  solid[2}  and  luid  [i-T]  has  .  .  p  nx 

k|<  central  stage  in  computational  continuum  mechanics  u(r,z)-u(r,*)4^r.  (3) 

■march  and  many  thoughtful  papers  [1.9]  regularly 

Wjcw  on  bis  subject  in  the  open  fiterature.  Heverthe-  Subsequently  we  w9  use  aa  overbar  to  denote  values  at 


■to  be  author  dares  believe  that  no  significant  im- 
imotflKnt  has  yet  been  offered  ever  his  original  toch- 
■We,  which  he  tecs  as  the  most  practicaly  satisfying 

*  21  *eoretieafiy  safe  method  to  cemputationeiy  deal 
**h«efid  aeariucompressMlity. 

5  Jh  technique  it  not  a  penalty  (the  abhor  hat 
[  *"  pot  reconcflod  himoctf  to  bis  atfy  tom) 
22*19-12],  w  im>  R  raquirt  reduced  (not  la  be 
25*md  with  numerical)  toasgration  h  it  a  standard 

*  ttiun  ef  the  displacement  finite  ohawnt  method  n 
.  **Ml  atomicity  bdudbg  a  wise  Imit  on  No  hub 
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v  ■  1/2,  as  •  ■  a(r,«). 

A  redial  displacement  a  i 
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K  being  the  buft  modulus  which  it  equal  to  x  +  2/3.  Abo, 
toe  stresses 

*r(r,  z)-t(l-pj»mi  «Ar,  (0 

are  both  independent  of  x. 

Equation  (3)  can  be  written  as 

mr,x)ma(r,m)*pT^-^^pT..^  (5) 


a(r,  z)-  a(r,  •)+()(*"  *)  (10) 

in  agreeacnt  with  the  general  prediction  of  Ref.[l], 

lUNutnmiMni 
Owr  sphere  has  stored  clastic  energy 

E"J [  (jaf,+et*4a,*+«*V^  (II) 

and  throughout  this  paper  wt  shall  assume  that  a  ■  1ft. 
For  a  single  element  that  extends  betweea  (»-|  and 
|>l(err«riandr«r,4h) 

£  -  f /’  [j  *(«< + 2»)*+(f*d*+2a^J  H  (12) 

where  dot  denotes  dUfcmndstous  with  resoecl  to  L  and 
where  g -M*. 

An  exact  stifness  antrix  for  a  two  aodal  point  afcacat 
is  reader  derived  (ran  eqa  (12)  in  the  (one 


i  iii  i  i  ■  ■rTufnp 


-  (IS) 

in  which  n  •  Ul*  rjk  ad  $  •  W+hlk.  a  ♦  0  •  1 


To  mderstaad  in  what  tray  this  element  lochs  whe, 
*-»•  observe  tost  that  linear  shape  function  product  t 

«- j(u,-ai)+ j(«,r,-ivO  (|g  , 

j 

so  that  «"0  is  possfek  only  with  Uir,-u»fi»l  am  i 
a]-ai>0,  or  u, »  a* •  0.0a  toe  other  hand  die  antsm  * 
error  in  toe  tonic  element  computed  A  A  and  to  is 

Btu-D-jsjfie-DVdr 

+ £  Ke,  -  tof  4  (e.  -  toflr*  dr  CtSA 

As  x-*»,  the  exact  r-*0,  and  because  o I  toe  large  x  o 
toe  dOatoric  energy  expression,  I  is  forced  down  too  sal 
locking  ensues.  Choosing  toe  tiro  parameter  iaitt- 

“*lc,r+3  0* 


for  the  element  displacements  allows  a  constant  dQata 
tion  without  locking  but  we  forgo  this  restricted  pot 
ability  here. 

We  now  emplo>  the  linear  dement  to  compute  tot 
elastic  displacements  of  the  radially  deforming  sphere  tf 
Section  1.  Since  we  use  the  exact  variational  principle  m 
are  assured  an  Ofli*)  convergence  in  the  displscaneau 
for  any  chosen  value  of  x.  because  of  the  factor  x  a 
£(«  -  d)  in  eqn  (15)  we  must  also  be  resigned  to  a  factor 
x  in  toe  displacements  error,  and  a  locking  mechanism 
Figure  I  shows  the  error  in  a,  ■  u(o  »  1/2,  x)  for  diffcrcu 
values  of  z,  and  indeed  for  a  number  of  elements  Nr 
suCdently  large  |tort/u,|  is  proportional  to  Ne'*.  Figure  2 
shouts  the  error  growth  in  a,  ns  i  is  increased,  in 
different  values  of  Nc,  and  truly  when  x  <  Ne*  toe  com 
puted  error  in,  toe  x  the  dependent,  a,  is  proportional  Is 
x.  Precisely 

f  ‘ 

|^|-0J2xNe'*  (P) 

and  we  conclude  that  for  any  fixed  value  of  x  any  desks!  | 


ghplacemeat  discretization  accuracy  cm  be  teemed  wife 
t  idficieotly  large  auaber  of  elements,  but  the  larger  z 
At  MR  elements  seeded  for  any  gives  accuncy.  A 
Ml  e ft- 1000  ces  be  safely  assumed  for  the  avaflabk 
tlattk  materials.  At  tbit  value  of  z,  r»  0.9995,  asd 
accordiag  to  eqn  (3) 

a(|.  10*)  •  u(j.-) +130 10-’.  (10) 

Eves  if  tbe  elastic  material  is  practically  incom- 
pressMe  it  would  be  unwise  to  substitute  a  targe  z  value 
is  the  faite  element  stiffness  matrix.  It  is  computation' 
sly  more  sensible  to  equate  the  relative  compressMity 
era,  which  for  ut  is  about  l/6z,  with  the  discretization 
enor,  for  a  smaller,  mesh  dependent  z.  Our  linear  de- 
meat  invites  the  choice 

0.32zNe‘,«0.!67z-1  (19) 

At  leads  to 

X  “  0.72  Ne  (20) 

asd  a  balanced  error  of 

|^|-0.23Ne'’  (21) 

•here  u,  is  a  function  of  z  and  consequently  Ne. 

At  this  point  many  a  reader  will  raise  the  objection 
that  the  optimal  z/Ne  ratio  is  not  only  element  dependent 
hut  also  problem  dependent  and  is  possible  here  only 
with  the  aid  of  the  exact  solution.  This  is  true  but  we 
mast  bear  in  mind  that  the  computational  analyst  is 
■eguhrly  confronted  with  the  question  of  choosing  dis¬ 
cretization  parameters  without  ever  being  able  to  select 
Ac  optimal  ones.  Does  he  ever  know  what  number  of 
dements  or  which  mesh  layout  is  optimal?  All  we  can 
do,  for  a  given  element,  is  settle  for  a  reasonable  z/Ne 
"tio  and  either  incorporate  it  in  a  Unite  element  com¬ 
puter  code  or  instruct  the  user  as  to  the  proper  limit  os  z 
depending  on  the  number  of  elements.  What  is  important 
her*  is  that  a  constant  z/Ne  ratio  assures  a  simultaneous, 
*•*•,  h-»0,  energy  convergence  OfNe-*)  to  the  incom¬ 
pressible  solution.  One  order  iff  Ne-1  is  lost  in  this 
process  but  we  are  for  that  on  theoretically  safe  ground. 

When  the  stiffness  matrix  for  die  linear  element  is 
created  with  a  one  point  Gauss  integration  rule  one 
dobnl  lute  clement  made  escapes  locking  and  the  dda- 
Anc  part  of  the  global  stiffness  matrix  is  reduced  from 
■dug  positive  definite  to  being  positive  lemidefinite  with 
***  mo  eigenvalue.  It  happens  for  the  sphere  that  this 
fdt  is  just  the  right  one  for  the  approximation  of  the 
■cumpmaMe  solution  and  finite  elements  produce  now 
cccoatc  nodal  values  independent  of  z.  This  is  the  trick 
*  "dated  integration;  to  save  (good)  approximation 
y°des  from  getting  locked.  The  success  of  this  con- 
**Aeut  device  appears  however  to  be  only  conditional 
y  »  very  large  z  can  leave  at  our  disposition  lanky 
JProximatiou  modes.  Counting  degrees  of  freedom  may 
the  number  of  locking  tm  modes  in  the  db- 
rAoation  hut  foretefls  nothing  about  their  quality.  Here 
"A*  weak  point  d  penalty  methods:  they  insist  on 
^1  taq  computatiooaly  pure  incomprescAlc  modes. 

i  —an— asmuini 

■  conclude  from  Ae  previous  section  that  for  aearty 
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order  of  the  element  for  better  computational  efficiency. 
When  exact  quadratic  elements  are  used  to  discretize  the 
elastic  sphere  the  relative  error  in  a(  is,  according  to 
Figs.  3  and  4 

|^|  -  0.034  4zNe“  (22) 

provided  that  z<Ne4.  The  compressibility  error  stays 
the  same  1/6  z,  and  therefore  incompressibility  is  bat 
approached  with 

z  -  2.2  Ne2  (23) 

rouhingtn 

|^|  -  0.076  Ne  *  (20 

where  n(  is  stfil  a  function  of  z  a  Ne. 
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Fn  S.  Error  in  die  ncooprettMe  i>  H  depending  oa  z  for 
different  number  of  finite  dementi  Ne  in  die  ditcrttization 


Fijure  $  clearly  shows  the  convergence  procesi  to  the 
totally  incompressible  solution  t,.  For  a  fixed  number  of 
elements  as  z  is  increased  the  approximation  is  first 
improved  as  the  compressibility  error  declines.  An  exact 
value  happens  to  be  crossed,  but  further  increase  of  z 
brings  only  a  decline  in  accuracy  as  the  element  stiffens 
and  ultimately  locks.  Figure  5  shows  also  the  error 
reduction  in  fi,  for  the  choice  z-4  Ne1.  Notice  that  with 
Ne  -  20  we  can  go  here  as  high  as  z  -  MOO  or  *  -  0.4997, 
and  get  a  relative  error  |88/8|  •  MT*4*. 


calling  for  round-off  precaution! [14,  IS].  Our  IBM  com¬ 
puter  hat  a  tingle  precision  relative  accuracy  of  Oi  MT* 
and  a  double  precision  relative  accuracy  of  I0~w.  On  this 
machine,  when  using  tingle  precision,  aO  accuracy  is  lost, 
as  toon  as  the  number  of  etc  menu  reaches  17.  Double 
precision  saves  the  day  for  nearly  incompressible  solid 
computations  by  banishing  this  unhappy  episode  to  Ne  » 
4500!  As  Fig.  6  shows  available  higher  computer  ac¬ 
curacy  permits  the  achievement  of  any  reasonable  dis¬ 
cretization  accuracy  practically  free  of  round-off  errors. 

«.  IXTluraLATIONS  TO  TRE  LMT 
Simultaneous  extrapolation  to  the  limits  of  z  and  Nt 
can  desist  the  need  for  their  large  values.  Let  8  denote 
Ae  incompressible  solution  and  a  the  Unite  element 
computed  displacement.  Since  the  compressibility  and 
discretization  errors  are  proportional  to  z~'  and  z,  res¬ 
pectively,  we  may  write 

a-ff+|+ffz  (27) 

and  three  a  values  for  three  z  values  will  determine  a,  & 
and  8.  Suppose  that  ui,  Uj  and  n>  are  the  three  a  values 
computed  for  z(,  z>  ■  2zi  and  z>  -  4z,.  then 

f--2ui  +  5ua-2a,  (2» 

ffii»|(ai-3a,+2aj  (29) 


We  infer  from  the  elastic  energy  expression  (11)  that 
the  global  stiffness  matrix  is  with  a  spectra)  condition 
number  [13]  proportional  to  z  and  Ne1.  Direct  com¬ 
putations  reveal  that  the  quadratic  element  leads  to 

Csm-fcbzNe’  (25) 

or  withz»4Ne*to 

Cj(K)  -  26.4  Ne4  (26) 


«/z,*  j(2n,-3uj  +  Uj).  (30) 

For  example,  the  quadratic  dement  of  Section  4  com¬ 
potes 

a, -0.14234053  for  z.  -  512 
■i-0.)4176086  for  Zj-1024  (31) 

■>-0.14064927  for  z,-2048 

with  which  eqn  (28)  yields  8-0.1428247  as  compared 
with  the  exact  8*0.1428571,  or  a  relative  error  of 
10-*4*.  Knowing  ffzi  in  eqn  (26)  allows  for  a  better 
approximation  for  the  compressible  a  at  z( 

a-j(2a,+3aJ-2aJ.  (32) 

Using  the  data  ia  eqn  (28)  we  compute  with  aqa  (32)  s 
■(1/2, 512)  -  0. 14288837  as  compared  with  the  exact  a(U2. 
512) -0. 14230359,  a  relative  error  of  Hr1*’,  the  con 
responding  a, -a(!/2, 512) -0.14234053  from  eqa  (28)  is 
with  a  relative  error  of  KT*4. 


We  forsee  dffkdties  in  the  computation  of  the  pres¬ 
sure  p  -  Ke  as  sec  have  to  multiply  a  smal  e  by  a  Imp 
K,  and  stB  expect  aa  accurate  p.  Fifarat  7  and  •  dner 
the  computed  pressure  p-Kc  over  e  aweh  of  32  bar 
laile  dements.  The  change  of  aiga  ia  p  that  takae  pieoe 
inside  each  dement  can  be  explained  whh  refarsnot  is 
the  energy  error  ia  eqn  (15)  and  its  mMutaadeu  by  tit 
Ante  dement  eolation.  We  expect  alia  a  paad  evemt  p 
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I  Hg.7.  Computed  faiun  dhtribuboe  m  cud  1 


r  points  at  the  optimal  integration  node*  for  a  weight 
j  Attribution  r*. 

For  the  linear  element  optimal  integration  point  it  at 
-  Ac  root  of  the  equation  r-b»~0  where  minimixes 


fc-  r(h.-r)r*dr. 


*  A(-j(r1,-r,1)  (34) 

.  we  have  that  in  the  interval  -1  £  f  s  1  this  point  it  at 


EnociMopma/  NJ\V\\ 

Urmcr  tttmtnts 

Ono  point  ktograOom''  'O'X 

Rg.  9.  Caevergcac*  of  avenge  preamn  m  the  hrst  dement 
(between  r,  - 1/2  and  n“  1/24  k)  with  Nt,  m  dcpeadag  on  x. 
Cmvesufor  exact  iategratioe  of  the  drviMoric  pen,  rod  one  pomt 
~*  rut  nf  tkr  total p niiatid  earngj  Pm  m  I  ■ 


HGH-' 


■mch  means  that  in  a  32  element  mesh,  in  the  flrst 
dement  ( - 0.01025562.  while  in  the  last  clement  (» 
M0524924.  Figure  9  describes  the  convergence  of  the 


_  '••£/>**  «• 

a  Ae  first  element  for  both  the  exact  and  approximate 


dement  To  locate  the  optimal  integration  points  for  this 
dement  we  minimize 

Ltmf  (he+hir-r^Vdr  (37) 

Jn 

with  respect  to  ft.  and  h.  and  Ind  the  roots  of  the 
equation  h«+hir-r*»0.  Doing  tins  we  have  for  the 
flnt  dement  between  r, ■  1/2  and  r*»  1/2+ 1/16, 

-  0.561056  and  fa  -  0.59247*  (31) 


^Qmdratic  elements  are  more  interesting.  Figure  10 
yw  the  computed  pressure  Attribution  over  an  I 
•**nt  mesh  of  quadratic  dements.  Change  of  sign  in 
■*  pressure  error  distribution  occurs  now  twice  in  such 


m 
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Tads  1.  CompcUd  average  prtuart  p',  and  the  pm  Hires  at  the  Gauss  poents  Gi  and  G:.  at  well  as  at  the  aptr-wl 
peats  Fi  and  Ft,  for  aa  t  element  discretization  with  2  »  254.  Exact  p  ■  0.MM57U 

elenset  p  a(s,)  »(f,)  MFj) 


1 

0.144)4 

0.1S3S4 

0.10137 

0.14)41 

0.1434* 

z 

0.14)40 

0.14044 

0.12120 

0.142*0 

0.14301 

3 

0.14ZS0 

0.14444 

0.13074 

0.1427* 

0.142*4 

4 

0.14280 

0.14048 

0.13446 

0.14270 

0.14272 

S 

0. 14271 

0.14747 

0.13037 

0.14244 

0.14244 

s 

O. 147(4 

0.14444 

0.13**2 

0.1424) 

0.14243 

7 

0. 14243 

0.14447 

0.14004* 

0. 14241 

0.14242 

e 

0.14242 

0.143*1 

0.34141 

0.14241 

0.14241 

at  compand  with  the  two  Gains  points  (,  j  -  z  0.577350. 
Table  1  lists  the  computed  pressures  at  the  two  Gauss 
points  G(  and  Gi,  at  the  optimal  points  Ft  and  f»,  and 
also  the  average  pressure  inside  each  ekatent  for  the 
mesh  and  z  d  Fig.  10.  Figure  II  shows  the  convergence 
with  Ne  of  Ac  average  pressure  m  the  irst  element  of 


Fin.  10. 

Figure  12  shows  the  computed  pressure  distribution  in 
a  four  cnbk  element  discretization  of  the 


when  At  is  given  in  oqn  (34).  bt  the  flrst  cable  dement 
of  fig.  12  dm  three  optimal  points  an  at 


(,•-0.76076446 

fc- 103 173007  (40) 

(,-0.7N37lll 

whanos  dteOenm  points  an  of  fij- a  0.77490667, 6" 
0.  At  the  opirnnl  intagnaioa  points  we  Ini,  lor  >  •  NM, 


the  three  corresponding  computed  pressure  values  p* 
0.129513,  0.161753,  0.127681.  At  the  three  Gauss  porno 
we  compute  the  less  accurate  p  ■  0.237300,  0.053M1, 
0.199470,  that  we  compare  with  the  exact  p  ■»  0.142857. 

» 

smavKon  t 

Our  present  discussion  wfll  be  wantiag  without  t 
mukidimenstooel  example.  Let's  considfr  for  that  pm- 
pose  the  problem  da  hollow  sphen  revolving  around  m 
axis  that  passes  through  its  center!  16).  Both  the  hnar 
radius  r  ■  1/2,  and  the  outer  ndins  r  - 1  d  the  sphsn  nee 
assumed  free  d  surface  (dees.  Symmetry  permits  is 
reduction  d  the  problem  to  the  portion  shows  fat  Fig.  & 
Discretization  d  the  sphere  is  sccompBthwl  with  6>t 
biquadratic  elements  nwncricdfy  hnegromd.  Aho  edh 
reference  to  Fig.  13  we  sesame  •  •Otobeteaxktf 
revolution,  and  suppose  for  this  exempts  x*23t. 

Table  2  lists  the  computed  avenge  present 
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Table  2.  Pressure  «l  the  Saite  etancats  of  Fig.  11  Vetoes  arc.  tram  lop  to  bottom  is  Mck  square:  excel  preswre  at 
Ike  center  of  ike  element,  avenge  over  clement  with  a  J  *  )  aSegntion  rale.  and  avenge  ever  element  with  a  1*2 
totearation  rale.  Square  at  bottom  left  ia  for  rime  I,  wkie  at  top  rgbt  it  it  for  eternal  36 


0.403300 

0.3(0233 

0.311(03 


0.370445 

0.360236 

0.360104 


0.256531 

0.247030 

0.24M30 

0.133455 

0.127514 

0.120054 

0.245(00 

0.245330 

0.23(t0t 

0.130760 

0. 13(601 

0.134175 

0.225775 

0.225071 

0.221424 

0.147(42 

0.140146 

0.143(01 

0.204700 

0.203030 

0.20007 

0.150545 

0.157(02 

0.15S26( 

0.1(5500 

0.104464 

0.11)072 

0.1(7727 

0.165374 

0.1(5031 

0.174500 

0.173(17 

0.172000 

0.373020 

0.171271 

0.170675 

0.01(543 
0 .017415 
0.016352 


0.040(11 

0.035537 


0. 070(50 
0.071(35 
0.060147 


0.1005(6 

0.100157 

0.105756 


0.1411(1 

0.13(4(6 

0.130316 


0.160010 

0.157507 

0.157111 


-0.007360 

-0.00556 

-0.00(103 


L-  f  (ot+Oir-rYr* triad 

UmL  <*•♦*•'*  rV,hl<# 


ft— 13SS7K.6-ISN906 


■  die  element  that  Set  between  r,  - 1/2  aad  r,  - 1/2 + 1/12. 


i»,- -0.2906563,  0.6192173  (44) 

ia  the  ome  element  that «  between  I, -  0  aad  #,  -  nf  12. 

At  these  four  integration  points,  numbered  si  is  Fig. 
13.  wo  compute  p  -  0.6185,  0.7719,  0.4360,  0.7439,  os 
compared  with  dm  cuct  p  -  0.7000,  0.7123.  04300, 
0.6333.  Computation  being  done  with  3x3  Gaum  ra¬ 
ted  elements  and  x  -  236. 

As  for  the  displacements  we  have  at  point  A 


-0.1009331  (z-256) 

-0.1007301  (i  —  (a) 

-0.1131709  (2x2)error-3J7%  '  ' 

-0.1069602  (3  x  3) error-  1J2%. 


Wa  may  attempt  to  approximate  the  pressure  at  point 
C  through  a  bffiaial  extrapolatiou  over  elements  23. 26, 
31  and  32.  Let  p»,  etc.  denote  the  average  pressure  at 
element  26  assigned  to  its  center.  Extrapolatiou  bora  the 
four  center  values  furnishes  the  pressure  approximation 
for  point  C 


P  “  j  (”3pm  ♦  p»4  9p»i  -  3pm) 


Pmm,- 0.72479  (3xJ) 

Pm. -0.74613  (a-23« 
9^-0.74067  (a— •) 

t  a  dberetiaation  error  af  about  3*. 
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1  I nt roduct ion 

Unconditionally  stable1  (semi)  explicit  integration  schemes  for  stiff  svstems 
of  equations  as  generated  by  finite  element  analysis  of  elastodynamics  and  non- 
stationarv  heat  transfer  are  shown  to  suffer  from  the  DuFort-Frankle-Saul ' ev  [l] 
syndrome  wtiereby  coupling  between  the  space  and  time  discretizations  may  have  a 
ruinous  effect  on  the  accuracy  of  the  computations. 


2.  Elastodvnamics 


To  integrate  the  discrete  equation  of  motion 


+K^j  =0 


in  which  M  is  the  global  mass  matrix  and  K  the  stiffness  matrix,  a  (semi)  ex¬ 
plicit  method  has  been  introduced  [2]  in  recent  years  into  the  engineering  liter¬ 
ature  (see  also  Refs.  ['1,4])  based  on  the  s\iui:fct  ric  splitting,  of  K  into 


K-L+L 


where  I,  is  lower  triangular.  Then  for  time  step  X 

(M+ 


(m+£lt)3i  =(m- 


:2,t 


“2- 


It  M  is  diagonal  passant'  I  rom  time  level  0  to  time  level  1  requires  back 
substitutions  only  and  lienee  the  appeal  of  scheme  (3). 

Tt  is  shown  [5]  that  this  method  is  energv  conserving  and  unconditionally 
stable.  It  is  also  shown  in  Ref.  [5]  that  for  any  fixed  !,  spacial  mesh  re¬ 
finement,  done  with  the  hope  of  improving  the  discretization  accuracy,  will  even¬ 
tually  sabotage  the  quality  of  the  displacements  computed  with  eqns  (3)  or  (4). 

We  reconsider  here  the  singular  behavior  of  the  computational  errors. 


! 


}• 


Consider  the  string  problem 


,U 


(jCt  o)  -  S/*  Zlift 


(5) 


for  which  we  know  that  u(x,t)  =  sin2itxsin2‘,'t .  According  to  the  analysis  of 
Ref.  [4]  when  eqn  (3)  is  discretized  with  first  order  elements  of  size  h  and 
a  lumped  mass  matrix,  the  first  mode  excitation  causes  the  string  to  vibrate  with 
the  computed  frequency  given  bv 


— - — — ^ 

4  +  3l  (i-cos2nt,)f  l(s+  'f'') 


(6) 


where  it  -  i/li.  Power  KcrlfH  expansion  yields 
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Co  (3  ^  +  H 


(7) 


(8) 


This  error  estimate  for  the  computed  frequency  includes  the  expected  0(h  )  from 

the  spacial  discretization,  O(t^)  from  the  time  discretization,  and  also  the 
A  2 

coupled  0([  /h  ).  Holding  (  constant  and  reducing  h  lowers  i5u>  until  a 

-k  it 

critical  vt'  =  (l/h)  is  reached,  beyond  which  lc  i nc reases  and  accuracy  is 

lost.  This  behavior  of  ooj  is  clearly  seen  in  Tig.  1  directlv  computed  from 
eqn  (6). 

The  critical  V  is  found  from  dw  /dh  =  0,  and  we  find  that 


i  h 


■0 


(9) 


One  would  not  want  to  exceed  this  sl<,  even  for  the  stable  method.  Going  beyond 
does  not  produce  the  computational  hang  typical  of  instability  but  its  effects 
could  be  as  insidious. 


Figure  2  shows  the  computed  motion  of  the  central  point  of  a  string  discre¬ 
tized  with  linear  finite  elements.  With  h  =  0.05  and  ^  =  0.1  the  computed  so¬ 
lution  can  barely  be  distinguished  from  the  exact,  but  small  elements  cause  a  hefty 


error  in  the  period. 
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3.  Heat  Transfer 


For  the  solution  of  the  system  of  first  order  equations 


the  unconditionally  stable,  explicit  scheme 

=  c  f  (*) 

5  -  3*. 

aA  -  M  ±  ($<  3) 

0!  <7°  CjT*j 


has  been  advanced  [7,8]  recently. 


To  analyze  the  scheme  in  eqn  (11)  we  consider  the  linear  system 


3+K3"° 


with  positive  definite  and  symmetric.  K.  Let  v^  and  v,;  be  two  eigenvectors 
of  K  and  0  «"  )  **  the  corresponding  eigenvalues.  The  initial  value 


V  ^ 


produces  the  evolution 


-M(f)  -dtt)^ 

+•  by) 

-  e 

b(0  «  e 

(14) 

Or  with  t  =  j'l,  j  =  0,1,2,... 

,  C'j  =  and  02  =  /.,i 

-j'H 

a-.  -  a. 6 

J  2 

J,-  =  beJ^ 

J 

(lb) 

Starting  with  y  in  eqn  (13)  we  have  that 

%  ~  -  ( T?,  +  b  ) 

and  eqn  (11)  produces 


Figures  3,4,5  and  6  describe  the  behavior  of  eqn  (17)» 

Starting  with  a  =  b  =  1  in  eqn  (13),  Fig.  3  shows  the  error  6a  in  a(t) 

* 

at  t  =  1  as  a  function  of  ^ •  It  is  clearly  seen  in  Fig.  3  that  ^  ~  2  con¬ 
stitutes  a  critical  value;  below  it  the  error  in  a(l)  is  nearly  independent  of 


4*2  but  as  ^2  *  ^  is  crossed  the  accuracy  of  a(t)  suffers  a  sudden  loss.  A 
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smaller  and  more  steps  to  reach  t  =  1  produces  a  more  accurate  solution 

for  <  2,  as  can  be  seen  in  Fin.  **  but  as  soon  as  =  2  is  crossed  the 

same  sudden  loss  of  accuracy  reoccurs. 

In  Figs.  5  and  b  we  examine  the  errors  in  a(t  =  1)  and  b(t  =1)  as  a 
function  of  the  initial  b,  when  the  initial  a  =  1.  It  is  readily  concluded 
from  Figs.  5  and  6  that  high  machine  accuracy  is  of  little  help  in  avoiding  the 
sudden  error  jump  in  a.  Only  when  b  =  0,  is 


(18) 
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Abstract 


An  alternative  to  the  Riks-Wempner-Crisfield  iterative  correction  scheme 
is  presented  that  does  not  require  an  explicit  displacement-load  accession  path 
to  the  nonlinear  equilibrium  curve,  nor  a  known  equilibrium  point.  Its  syametry 
with  respect  to  the  displacement  and  load  assures  success  in  rounding  limit 
points  as  well  as  turning  points. 
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1.  Introduction 


Consider  the  single, implicit  equilibrium  curve 


rc#,?i)  —o 


m 


in  which  x  is  the  displacement  and  X  the  load.  To  trace  X  versus  x  we 
shall  have  to  compute  close  root  pairs  (x,  X)  to  eqn  (1). 

Since  Iterative  methods  must  Invariably  be  used  to  solve  the  nonlinear 

r(x,  X  )  »  0,  two  strategy  decisions  concerning  the  solution  procedure  have  to  be 

made:  (i)  how  to  advance  from  an  established  equilibrium  point,  say  (x  ,  X  ),  to 

o  o 

a  next  Initial  say  (x^,  X^);  and  (11)  what  corrective  method  to  use  if 

r(x,s  Xj)  is  untolerably  large. 

The  simplest  and  cheapest  tracking,  or  continuation,  procedure  for  r(x,  X)  -  0 
consists  of  increasing  Xq  to  Xj  for  a  new  Initial  guess  -  xq,  X^)  and  a 

(modified)  Nevton-Raphson  Iterative  solution  of  r(x,  Xj),  with  a  constant  load 


That  a  monotonic  X  sequence  can  miss  sections  of  the  equilibrium  curve  has 
long  been  pointed  out  [l,2].  To  remedy  this  Rika  [3,5],  Wempner  [4],  Crlsfleld  [6,7] 
and  others  [8,9]  suggest  an  iterative  approach  to  the  equilibrium  curve  on  a  var¬ 
iable  X  ■  X(x)  load-displacement  constraint.  The  difference  between  the  different 
methods  being  that  while  like  and  Vempner  advocate  a  linear  (planar)  constraint, 
Crlsfleld  promotes  a  spherical  one. 


Continuation  [lO-ll]  of  the  equilibrium  curve  is  invariably  baaed  on  linear¬ 


isation.  Let  (xQ,  Xfl)  be  a  point  in  the  (x,X)  plane  not  necessarily  satisfying 

the  equilibrium  equation.  r(xn,Xn>  i  0.  Let  further  fix  •  xn) ^  -  xq  and 

6X  ■  X  -  X  denote  the  displacement  and  load  corrections.  Linearisation  of 
0+1  n 

r(x  +  fix,  X  +  fi,)  ■  0  yields 
n  n  a 

-o  « 

where  r'  “  dr/dx  and  r  m  dr/dX.  Henceforth  we  shall  omit  the  subscript  n  when 
referring  to  the  current  nth  point. 

As  for  the  prediction,  that  is  for  the  move  from  the  already  computed  equi¬ 
librium  point  (xo,  Xq)  to  the  next  initial  guess  (Xj,  Xj),  it  is  commonly 
agreed  to  supplement  eqn.  (2)  with  the  elliptical  constraint 

o?  <^a~sx  <3> 

where  a  and  6  are  scaling  parameters  and  s  the  step  size.  Combining  eqns.  (2) 
and  (3)  we  get 

|  » 

since  r(x  ,  X  )  •  0.  The  choice  of  sign  In  eqn.  (fi)  determines  the  direction  of 
o  o 

travel. 

In  this  note  we  present  and  experiment  with  a  correction  procedure^ based  on 
the  Hewton-Iaphson  method  a a  axpreased  In  aqn.  (2),  whereby  the  initial  guess 
(Xj,  Xj)  la  iteratively  improved  on  an  orthogonal  trajectory  to  the  equilibrium 
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curve,  but  with  no  analytical  tortMlon  for  It  and  without  the  Involvement  of  the 
previously  computed  Xe). 

2.  The  Rlka-Weapner  method 

The  cleverness  of  this  method  lie*  In  the  replacement  of  the  constant  X  cor¬ 
rection  by  the  linear  relationship 

*  %  *  j C  (5) 


so  that  a6x  +  $6X  ■  0.  Adding  this  constraint  to  eqn.  (2)  yields 


pc  r 

Pr'-ocr 


pr 

P  r'-ec  r 


(6) 


or  If  the  constraint  line  in  eqn  (5)  Is  chosen  to  be  orthogonal  to  the  tangent  vec¬ 
tor  (6xo,  6Xq) ,  as  Rlks  and  Weapner  suggest,  then 


r 


(7) 


In  preparation  for  things  to  com  we  further  observe  thst  since  r^  6xq  +  rQ  6Xq 
eqn.  (7)  My  be  rewritten  as 


rr0 


(8) 


fYJ  +irfi 

If  r(x ,  X)  ■  f(x)  -  X,  then  r'  ■  f* ,  r  ■  -1,  and  eqn.  (8)  bec< 


0. 


(9) 
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that  is  shown  to  be  directly  obtainable  froa  the  Newton-Raphaoh  method  applied  to 
the  intersection  of  r  -  r(x,X)  and  X  -  -x/f^  +  &.  Indeed ^slnce  fix  «-r/r* ,  we 
have  that  fix  •  -r/(f*  -  ),  but  X*  ■  -Vf*  and  hence  eqn  (9). 

For  the  success  and  failure  of  the  Rika-VJ  emptier  method  in  turning  a  limit  point 
see  Fig.  1. 

3.  The  Crisfield  method 

Crlsfleld  replaces  the  linear  constraint  in  eqn  (5)  by  the  circular  (ellipti¬ 
cal) 


Using  the  linearization  r  +  r'  fix  +  rfiX  -  0  and  the  constraint 

(x  +  fix  -  x  +  (X  +  6X  -  X  5X  is  eliminated  and  we  are  left  with  the 

o  o 

quadratic  equation  [12] 

S  * V  f'\-  r)  {.«' +  CwO  f -  0  -a.)  f  f ' ] 

ail 

+  rx-2(?-»»)rf  -o 

for  fix.  The  lot,t{  increment  fiX  la  obtained  either  from  r  +  r  fix  +  rfiX  ■  0  or 
directly  from  the  circular  constraint. 

2  2  2 

A  direct  Newton-Raphaon  application  to  r  ■  r(x^X),  (x  -  xQ)  +  (X  -  Xo)  •  a 
produces  fix  ■  -r/dr/dx,  dr/dx  »  r*  -  rdX/dx,  (x  -  xQ)dx  +  (X  -  XQ)dX  •  0  OmJL 

h  - - - -  <12) 

fL  r  *-*• 

a -a. 

or  alternatively 


-’AiSRiTr- 
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r 

r  7 

It  is  interesting  to  compare  eqns.  (12),  (13)  with  eqn.  (7). 

For  the  success  and  failure  of  the  Crisfield  method  in  rounding  limit  and 
turning  points  see  Fig.  2. 


4.  Orthogonal  trajectory 


What  we  consider  now  is  an  orthogonal  trajectory  approach  to  the  equilibrium 
curve  without  the  need  for  an  explicit  constraint,  and  without  the  involvement  of 
an  equilibrium  point  (xq,  Xq).  All  we  have  to  do  for  that  is  apply  the  Newton- 
Raphson  method  to  r  -  r(x,  X)  -  0  so  as  to  have  6x  ■  -r/dr/dx,  add  to  this 
dr  *  x  dr  +  rdX,  and  Impose  the  orthogonality  condition 


(14) 


and  we  get 


_ 

r 


r*  +r 


•  X 


(15) 


which  is  eyvetrlc  in  x  and  X.  Comparing  this  with  eqn.  (8)  for  the  Rlks-Wemp- 
ner  method  reveals  how  the  orthogonal  constraint  is  updated  here. 


In  fact,  lean  [13]  has  suggested  a  correction  method  that  includes  updated 
normals.  Whan  applied  to  r(x,  X)  ■  0,  Rama's  method  produces  the  load  increment 


r 


06) 


involving  not  only  the  previous  equilibrium  point  (x  ,  X  )  but  also  the  new  ini- 

o  o 

tlal  guess  (x^,  X,).  It  is  interesting  to  compare  eqn.  (16)  with  gqns.  (7)  and 
(13). 

5.  The  vector  form 

Our  main  Interest  is  systems  of  nonlinear  stiffness  equation  as  produced  by 
finite  elements  applied  to  elastic  problems,  and  we  turn  our  attention  in  this  sec¬ 
tion  to  the  application  of  the  orthogonal  trajectory  method  to  problems  with  many 
degrees  of  freedom. 

Let  ir(x,X)  be  the  total  potential  energy  of  the  discretisation,  with  x 
being  the  displacement  vector  and  X  the  scalar  load.  Here  r(x,X)  **  3ir/3x  is 
the  gradient  vector  of  x.  If  (xq,  Xq)  is,  again,  an  equilibrium  point  so  that 
r(xQ,  Xq)  »  0,  then  a  linear  expansion  around  it  is  written  as 

!£(*-*.)+  =  0  (”> 

where  dr/9x  is  the  stiffness  matrix,  say  K,  and  where  dr/dX  is  the  load  vec¬ 
tor,  say  p.  For  how  to  compute  K  and  p  see  [14-19]. 

A  prediction  with  step  else  s  subject  to  the  constraint 

C*,  -%$(.%,-%.') + £»,-  2.)1-  s*  do 

together  with  the  linearisation 


(19) 
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leads  to  the  Initial  guess 


(20) 


where  6A  ■  A,  -  A  . 

oxo 

Newton-Raphson's  correction  applied  to  the  system  of  linear  equations 
r(x,A)  -  0  becomes 


tof) 


(21) 


.-1. 


A  constant  load  iteration  with  6A  ■  0  reduces  eqn.  (21)  to  6x  ■  -K  r.  A  lin¬ 
ear  load-displacement  constraint  that  relates  6A  to  6x  by 


£>  *  <£"  i* 


(22) 


with  the  given  vector  a  produces  the  correction 


ft- 


T  "I  y 

<lK  ^ 

i  +  <iVf 


,  t*— KV+tyO 


(23) 


For  orthogonal  trajectory  accession  we  K  •'f*  j 


(24) 


Observe  that  both  the  Rika-Wempner  method  and  the  present  one  call  for  the 
computation  of  K*lr  and  l”lp,  aad  otherwise  entail  comparable  computational 
effort. 


■*vmm***''*» 


6.  Midpoint  correction 

A  Runge-Kutta  type  correction  in  the  form  of  a  midpoint  reevaluation  of  the 
stiffness  matrix  and  load  vector  is  also  often  used  [20]  in  nonlinear  computation¬ 
al  elasticity.  In  this  method  an  initial  guess  is  predicted  along  the  tangent  to 
the  equilibrium  curve  with  eqn.  (20)  but  then  K  and  p  are  reformed  at 

x  »  x  +  6x/2  and  A  »  A  +  6X/2,  and  the  new  data  is  resubstituted  into  eqn.  (20). 
o  o 

Each  step  requires  in  this  way  two  assemblies  and  two  inversions. 

An  alternstlve  of  equal  computational  effort  would  be  a  linear  prediction 
with  one  Newton-Raphson  correction.  It  happens  that  this  procedure  is  considerably 
more  accurate  than  a  midpoint  correction. 

It  is  sufficient  that  we  analyze  the  methods  on  the  parabolic  stiffness  equa- 

2  2  2  2 
tion  A  *  ax  ,  where  a  “  A"/2,  with  the  circular  constraint  x  +  A  ■  s  . 

Their  exact  intersection  point  (xe,  A#)  is  at 

xt-k  [-£  M«V)T  ,  = £[-)  ♦(«*&?*]  <I5> 

Without  correction  point  p  in  Fig.  3,  (x  ■  s,  A  »  0^  constltutea  the  ap¬ 

proximation  to  equilibrium.  Since  both  the  exact  solution  (x#,  A#)  and  all  other 
approximations  are  at  an  equal  distance  s  from  the  origin  we  prefer  to  use  the 


directional  error 


if  *  cosT1  (26) 

-1 

of  (x,  A).  For  point  p  we  have  that  f  »  cos  (x*/s) ,  or  the  (expansions 

*V+J-*V+  — ) }  .••)  ti7) 
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ii 


get  that  for  prediction  only 

(f^*$ 


if  s  is  small. 


Midpoint  correction  reaches  point  m  of  Fig.  3  with 


(i+«V)4  > 


it 


0+*V)± 


Using  eqn.  (27) ,  and  with  similar  expansions  in  powers  of  s  of  x  and  X  in 
eqn.  (29)  we  have  for  small  s  that 

tf  -  £  *V  (30) 


which  is  a  substantial  improvement  over  eqn.  (28). 

One  Newton-Baphson  correction  with  a  circular  constraint  leads  to 


5  (14  3 «_<.)**  A=-XS+2-*sS* 

r+  (V  «4  1 


and  a  directional  error 


which  Is  two  orders  better  than  eqn.  (30). 


7.  numerical  examples 


The  ability  of  orthogonal  trajectory  accession  to  succeed  around  limit  points, 
bifurcation  polnta,  turning  points,  snap- through*  and  detachments  is  numerically 


*t**~-. 


demonstrated  on  two  simple  structural  elastic  problems. 


First  to  be  considered  Is  the  two  member  hinged  truss  shewn  In  Fig.  5  in  its 
original  undeformed  state.  As  the  load  X  Is  Increased  the  two  elastic  bars  shor¬ 
ten  symmetrically  until  a  snap-through  occurs  and  the  triangle  is  inverted. 

A  displacement  x  causes  the  strain 

e  «/-(•*+  %x-AvT  <33> 

In  the  bars.  With  each  bar  being  of  the  stiffness  k/2  the  total  potential  energy 

2 

of  the  two  bar  system  Is  ir(x)  ■  ke  -  Ax,  and  hence 

i,  r—4  ,  f'- kCe^+ee")  <*> 


where 

l-t 


I 


(35) 


The  exact  equilibrium  curve  r(x,A)  -  0  is  shown  in  Fig.  4. 

Figure  4  shows  also  the  convergence  paths  of  orthogonal  trajectory  accession 
given  by  eqn.  (15),  for  different  starting  points  that  are  far  from  equilibrium. 

To  observe  the  orthogonal  trajectory  one  must  look  close  to  the  end  of  the  iterative 
approach.  A  starting  point  on  the  normal  to  a  limit  point  converges  in  one  step. 
Notice  also  the  special  attraction  of  the  limit  points. 

Figure  5  shows  the  point  by  point  continuation  of  the  equilibrium  curve  of 
Fig.  4  with  the  predictor  of  eqn.  (4)  and  the  corrector  of  eqn.  (15),  with  3  New- 
ton-laphson  corrections  per  step. 
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The  two  member  truss  shown  in  Fig.  5  is  fist  when  X  ■  0.  Its  behavior  under 
load  Includes  bifurcation,  turning  points,  and  detachment.  Two  degrees  of  freedom, 
x  and  #,  determine  the  configuration  of  the  deformed  system,  whose  total  poten¬ 
tial  energy  reads  now 


- 1  cM  * 


where  kj  is  the  stiffness  of  each  bar,  k2  the  stiffness  of  the  torsional  spring 


at  the  vertex,  and  e  the  strain  in  the  bars: 


e  - 1  - 


/--y 


The  two  stiffness  equations  are  obtained  from  ir(x,l)  through  differentiation 

9 ft  —  U  fc,  4-C.  XT  *b 

gf  df  (38) 

n  * 

If  kj  -  1/2  and  k2  ■  e/16,  then  after  eliminating  x  between  the  pair  of 
eqns.  (38)  we  are  left  with 


When  w  ■  0  we  have  either  A  ■  x  or 


so  that  for  e  <  1  two  bifurcation  points  occur  on  the  X-axis,  and  they  coalesce 
for  c  ■  1.  As  e  surpasses  1  the  nonlinear  bifurcation  curve  detaches  from  the 
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X  axis  and  a  turning  point  (i.e.  dX/dx  «  °°)  la  created.  The  critical  <f*  and 


X  at  which  turning  occurs  are  given  by 


"Loim  <j> * 


+' 


=  6  - 


Jli  e>sj>* 


(41) 


For  a  given  e  >  1  no  nonlinear  solution  exists  for  whatever  X. 


Tangential  prediction  plus  orthogonal  trajectory  accession  has  no  difficulty 
,  following  the  equilibrium  curves  past  the  turning  points,  as  can  be  seen  in  Fig.  6. 

The  two  values  of  e  ■  1.0235  and  e  ■  1.1027  in  Fig.  6  correspond  to  <p*  ■  ir/12 
and  <fi*  “  tt/6,  respectively. 

Figure  7  uses  4  and  X  from  Fig.  6  to  trace  X  res.  x  using 

(jCelf  -/)  (42) 

Finally,  we  compare  the  performance  of  the  various  continuation  methods  dis¬ 
cussed  in  this  paper  in  rounding  a  limit  point  (dX/dx  -  0)  and  a  turning  point 
(dX/dx  “  °°) .  The  parabolic  r  ■  ax(l-x)  -  X  ■  0  has  a  limit  point  at  x  -  1/2, 

X  ■  a/4  and  is  of  the  shape  of  Fig.  1.  We  choose  a  ■  8,  xq  -  0.4375  (at  which 
|  r  ■  1),  and  X  ■  1.96875.  A  tangential  predictor  of  step  size  s  -  0.2  lands 

i  ° 

j  us  at  x.  »  0.57892136  and  X.  ■  2.1101714.  This  is  approximately  the  situation 

1  1 

I  shown  in  Fig.  1(b)  and  Fig.  2(a).  goth  the  linear  constraint  method  and 

i  Ream's  method  fall  to  find  an  equilibrium  point  starting  with  these  x^  and  X^. 

I  Only  the  circular  constraint  method  and  the  orthogonal  trajectory  method  converge 

here  -  to  different  points,  though,  as  can  be  seen  in  Table  1. 

Reversing  the  roles  of  X  and  x  :  r  ■  aX(l-X)-x  creates  an  equilibrium 


curve  as  in  Fig.  2(b)  with  a  turning  point.  Starting  from  the  reversed  Xj  *  0.57892136 

\ 

!  i 
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and  Xj  "  2.1101714,  the  circular  constraint  method  fails,  but  the  orthogonal  tra¬ 
jectory  method,  because  of  its  Inherent  symmetry  to  x  and  X,  converges  to  the 
reversed  values  of  x  and  X  in  Table  1. 


Circular  constraint  s  ■  0.2 

orthogonal  tra" 

jectory 

step 

X 

X 

X 

X 

0 

0.57892136 

2.1101714 

0.57892136 

2.1101714 

1 

0.62723680 

1.8891614 

m 

2.0485029 

2 

0.62287506 

1.8793660 

0.50023522 

2.0000049 

3. 

0.62283738 

1.8792878 

0.50023520 

1.9999996 

\m 

0.62283737 

1.8792878 

0.S00HS20 

l.wtllt 

Table  1:  Convergence  of  the  circular  constraint  and  orthogonal  trajectory 
methods  upon  r  ■  8x(l-x)  -  X. 
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Figures 


1)  Success  (a)  and  failure  (b)  of  the  linear  constraint  Method. 

2)  Success  (a)  and  failure  (b)  of  the  circular  constraint  Method. 

3)  Midpoint  and  Nevton-Raphson  corrections  with  circular  constraint. 

4)  Orthogonal  trajectory  convergence  to  the  equilibrium  curve  of  the  truss  in 
Fig.  5  from  different  starting  points. 

5)  Pointvise  tracking  of  the  equilibrium  curve  with  a  tangential  corrector  and 
orthogonal  trajectory  corrections. 

6)  Computed  equilibrium  points  for  the  two  member  truss. 
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Abstract 


A  cubic-cubic  element  stiffness  matrix  and  element  gradient  are  derived  for 
the  thin  shell  of  revolution  undergoing  large  axi symmetrical  Kirchhoff  deformation. 
Application  is  made  to  follow  the  nonlinear  elastic  distortion  of  a  spherical 
shell  under  surface  pressure  and  polar  forces. 
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1 . _ Introduction 

Discrete  integration  of  tiie  total  potential  energy  is  used  to  derive  the  stiff¬ 
ness  matrix  and  gradient  of  a  cubic-cubic  element  for  the  thin  shell  of  revolution 
that  may  undergo  large  elastic  deformation.  This  element  need  not  be  small  in  the 
sense  of  approximating  arc  length  and  is  routinely  assembled  for  a  Newton-Raphson 
solution.  The  same  technique  has  been  previously  applied  with  success  to  study  the 

large  deformation  of  the  circular  plate  [1  ],thc  cable  [2], the 
elastica  [;s],  the  rubber  membrane  [f],and  the  curved  extensible 
rod  It,  ]  . 

An  example  is  made  whereby  a  thin  sherical  shell  is  deformed 
under  the  combined  action  of  polar  forces  and  a  surface  pressure. 


2.  Elastic  energy 

Considering  axisymmetric  Kirchhoff  deformation,  shear  is  absent  and  the  elastic- 
energy  of  the  stressed  shell  reduces  to 

c  cc 

6  -  t  (07fe,+<iIfc2.-»-<rje3)ci'/  <0 

J  J  J 

V 

where  0  ,  are  the  normal  stresses,  and  ■  j,  >  the  corresponding 

strains.  By  v  we  denote  the  original,  unde  formed ,  volume  of  the  elastic  solid.  We 
further  assume  that  the  elastic  material  obeys  Hook's  law  so  that  stress  is  linearlv 


related  to  strain  by 


=  177 


¥ 

/ 
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where  K  is  the  elastic  modulus,  «*  Poisson's  ratio,  and  e  =  2  r  3  * 

lowing  the  common  assumption  that  o  =0  we  have  that 


(3) 


and 


(4) 


so  that 


£ 


i-p' 


(^1  **  ^2.  ) 


(5) 


in  case  of  a  rod  we  add  the  stress  assumption  =  0  and  have  that 


,  63=-^6, 


(6) 


and 


(2)  reduces  to  < t ^ 


1* 


The  clastic  energy  (5)  becomes 


(7) 


g 


fcEJ 


z,  E  \A(>)  0I3 


Fol- 
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wherc  s  denotes  arc  length  and  A  the  area  of  the  cross  section. 

3 .  Displacements  and  strains 

Let  (x,v,z)  be  a  point  on  the  middle  surface  of  the  deformed  shell  of 
revolution  so  that 

x  =  t-ty) 

2  2  2 

Obviously  x  +  y  =  r  .  It  is  helpful  to  introduce  the  angle  4>  measured  be¬ 
tween  r  and  the  tangent  to  the  generating  curve.  With  this  the  position  vector 

— ►  — ► 

r  and  unit  normal  to  the  middle  surface  n  are  written  as 

^Orws^rsl^z) }  n  =(sm<j>cos6>;jtn^s,»,^.-  <*s<£)  (,n) 

— ► 

The  position  vector  of  a  point  on  n  at  a  distance  c  from  the  middle  surface  is 

(11) 


and 

Jlik  -  otr  ■+■  ■+  £  ot/i  (12) 

■*+  — ►  — ►  -4  —4 

Since  dr  is  in  the  tangent  plane  and  n  •  n  =  1 ,  we  have  that  n  ■  dr  *  0, 

— >  2  h ► 

n  •  dn  =  0,  and  ds  =  dR  •  dR  becomes 


<1* *  olr*  4  dl?  •  d.i l 


(13) 
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where 


dr  ^  d»  +  i£  cte 

$2  &  9© 


otw  - 


2£  oL4>  + 

3f  r  9© 


9T_  _ 


*Z  9©  ' 

we  are  left  with 

o^r .  A\r  -aoc  A'£’+}cd9 


are  orthogonal 

-% 

— • 

f\  i- »  • 

>  W 

9© 

.  . 

9? 

>  9<|> 

9© 

2^ 

-a.  cj> 

f  2_  <j 

dLf.^/T  «  tt<j>C*g  + 


where  prime  denotes  differentiation  with  respect  to  t?  and 

t  t\“£  ^*1 

«-(*•' +2' J  (  Smcf*.—  ,  tos^.^  — 

<  zV-zV" 

*■  — — 


Final ly 


^ 4-') V + r s  1  n  ^ 6,1  +  ^ 
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as  compared  to  the  original 


c^o  ~  t*o +$()>')<£%+ (f0 


(19) 


Notice  that  no  distinction  is  made  between 


here 


(X 


«»  +  £«#J 


and  r  Hence  the  strains  are 


(20) 


and 

f-fo  +  — gitt4v)  (21) 

v; + jsiv)4>0 

4.  Integration  with  respect  to 

From  Eq.  (19)  we  get  the  element  of  volume 

=<i#.(i+sfr)(.u-$fl#)  Jxtyte  (22) 

in  which  =  <t>'  q/"q  and  =  sin<>  /ry  Substituting  dv  in  Eq.  (22)  and  the 
strains  K  j  ( O  and  in  Eqs.  (20)  and  (21)  into  the  elastic  energy  expres¬ 

sion  in  Eq.  (5)  we  integrate  it  analytically  with  respect  to  between  -t/2 
and  t/2,  where  t  denotes  the  thickness  of  the  shell. 

In  doing  this  we  encounter  integrals  of  the  form 
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4 

z 


\  h/Vf***  V  ^  /v3  ,  *n  -** _ 2-  (£p- 

J  *jrn  =  2^  *  +— ; 

-t 

•a  j- 


fll/l  4/1*f 


X 


,  <L(a-tlf-r)+w-m  0  t 

Jr  - - — - - - W,(AvXf|) 

<L 


b_ 

2, 

4 

2. 


(23) 


or  if  !axl  *  <■  \ 


Lfo  +t)  =  fl-x-f  <tV+ j  a?  x  -  f  <lV +0(«fi/) 


(24) 


and 


4l 


4 

2, 


/M  xV  /)  x%-  Px~4  V  __ 


a-X-M 


(25) 


Now 


g  =  £££.  j. 

^  .  M  #1 


•f 

where 


Z^fCK-x.]  1*^ 


/  J 


€i  -  **  ‘  >  r;  ’>*'"■  — — 

*o 

«,•  si±iit ,  ft .  iL , « .  iti 

/*  '  Of-  ’  ^ 


(27) 
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Tf  It  fi  /2!  and  It  r^/ 2 !  are  ignored  relative  to  1  the  elastic  energy  of  the 
shell  simplifies  into 


E- 


zxe  ± 

Z 


4 


(28) 


For  the  spherical  shell  Eq.  (28)  is  obtained  without  simplification  but  because 


:0  V 

5.  The  plate  and  beam 

When  the  shell  of  revolution  is  a  flat  circular  plate 


4o  ~  °  ;  Si  4*0  —  O  Z0  -  O  Z  p  r  o 

,  (29) 

*i- 2  >  r*~-{  >  *°  =  i 


and 


&=  frjz 

+  +  ^ +**'*' IT 

9  9  1 . 

where  r  =  (r*  +  z'  )^-l  and  »"  =  r/r()  -  1. 

For  the  beam  we  set  in  Eq.  (2b) 

^  s.1 


(30) 


(31) 


If  the  axial  extension  is  small  we  may  set  ■.  =  u  =  1,  and 


will) 


z.  z  X 

£, -(y+Y'1 )  -i 


6.  Stress  resultants 


,  *  = 


zV-  zVH 


(33) 


(34) 
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H,-  -^p — p  x-  £7.(^0- foi\ 


(37) 


7.  Differentiation  with  respect  to  a  vector 

Let  f (x)  =  f (Xj ,x^ , • . . ,xn)  be  a  scalar  function  of  the  vector  argument  x 
We  define  the  differentiation  of  f  with  respect  to  the  vector  x  as 

a£  _  r'_  «P  li  =  jp*_  1L_ 

W  *  9X,-  >  <9yx  *  n8) 


t  it 

Notice  that  f  i»  a  vector  and  f  is  a  symmetric  matrix. 
Obviously 

(f. f)'=^'  5  (ftf-tf'+fy' 

but 


(39) 


where 

pV\  2t  21 

+  <1  3X,  9*; 


is  a  nonsymmetric  matrix.  The  matrix 
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(42) 


is  symmetric. 

8.  Cubic-Cubic  shell  element 

Let  V  in  Eq.  (26)  measure  arc  along  the  generator  so  that  =  1  and 
d»?  =  ds.  The  finite  element  extends  between  Sj  and  Sj  +  h  so  that 
S’Sj+hS  0<?5l,  and  ds  =  hd? ,()',()  =h2()  .  To  have  a  cubic- 

cubic  C  element  we  choose  the  nodal  values  vector 

and  interpolate  r  and  z  with 


(4  3) 


///•I 


where  <t>  and  ^  are  shape  funefion  vectors 


(45) 


with 


(46) 


We  integrate  the  total  potential  energy  by  sampling  at  the  three  Gauss  points 


(47) 


-li~ 


with  weights 


w,-.w3-  | 


> 


The  values  of  r,z,r,z,r,z  at  the  three  integration  points  £.,  j  =  1,2,3 


are  given  by 


;  ^-£4,  , 

J 


where  8^  and  ^  shortly  stand  for  0(£.)  and  and 

</*'  I  -•  mW-SO  t  !.‘,V  »>.  5  •  \  ,vs.  0.  0,  50.,  1  >\  |S  ^  «.  (l  (|. 

<!>'.  k(4.  I,  o,  (i.  i.  |_  n  o) , 

•A:.'  '  !'il  (\  2  '  \  is.  I).  I).  0,  V  \  |S,  (I.  0) 

• h :  -  <\  1,0.  0,  0,  1. 1),  ()) 

-A',,  •  '  5  4sV  !>.<»,  0.  •  o\  15.  5?  3\  15.  0.0,. 

'A? ~  (<>,  l  o.  o.  o,  i,  o,  o) 


>0  ’  l.)V  15.  5  t-  V  I 

5,  0.  0.  50  i-  |J\  |S.  -  s 

■//, 

*(<>.  0.  4. 

1.0. 0.4.  1). 

*A'.» 

--  0. 

0.  ?  *  \  15.  t|,  0,  6. 

1  V  \1  5 ) 

'/»'■ 

'.(«».  0.  < 

'•  1.0,  0.6.  |). 

•-  1(0.  0.  . 

(A  15.  -  5  4  3V  Is. 

o,  0.  6\  |5.  5  3\  1  s\ 

'//'•  - 

(0.0,0. 

■1.  0.  (1,  0,  |) . 

of  h*>  is  for  Gauss  point  1  and  the  lower 


where  the  upper  sign  of  fl?  is  ror  bauss  puxm-  x  - - — 

for  Gauss  point  3. 

We  choose  to  derive  the  element  data  from  the  simpler  energy  expression  (28) 
and  write  for  t lie?  clli  element 


el-klZZ*;  Kj  (.^  .j  fe.j +  ) 

0*1 

+  (*lj  4-21/X^  W»j  +  *Wj  )] 


■■■■  “a?*  M  -  ■ 


ppiMP  «  AH  I  p*  If  III  Will 
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where  j  refers  to  the  jth  Gauss  sampling  point,  and  where 

.5.  UOi-r)  . 


b%  >  £ 


(52) 


From 


3e" 


9£e_ 

d^t 


»  * 

j?£e_ 


(53) 


we  have  that 

Jt 


$*•  =  i?,  ^  1I*J  [  ° L6 'J  S  *  ^S  €lJ  +  Si  e*j ) +  %  6,j  j 

+  *U  K-j  S  +  «!j  K,j  )  +K,J- 


(54) 


and 


*tm  \h  S  +  S' 4' **  ^S' V  S'  + 

*} .  x-,T  s  "  »  fT  M  ,  ,T  , 

S  S'  *SS‘H  +SS'+KD'S'+KijK*j+K»jS' 

*  >>  CKiiKlJ  +  Klj  Ky  ■*  Kij  K,J-  J-  Ky  Kj  )  J 


(55) 


where  (  )’  =  0/<hi 


To  shorten  the  notation  we  introduce 


f"  •  .  •»  .  ^  .  I, 

=  z:r-z.r  ,  +2 


(56) 


so  that 
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£2  =  nr- 1 


-i«-r  . 


.  -  r°  3  ^  < 


t'-=.yrnji  +<$•?■- ff-ztf  ,  ^  ~  ^(.r<f>  ■>■ 

f  ,,-$<^+4>fr~'j>$T-4>i''T ,  +yi'T) 


Next  we  write 


eXTg  lg'  ) 

£i  =Y'»',4  .  K,'--v--Vq-^/-  io-i 


which  is  all  we  need  to  program  the  element  gradient  and  matrix. 


9.  Potential  of  pressure  and  point  force 

i 

If  p  denotes  external  pressure  and  F  an  apex  ( r=  0)  force  then  tot 


we  add 


#  ,  * 
7C  -  1° 


rtli  -  f  z 


^  -A 'J'." 


* - ^4 

■  .  jg?  ■ 


and  have  that 


3.  =  f*  t  SwiC2^jt+vi^) 

j  -i 


(63) 


and 


k  =  z4-  2  wi  L£j  ■ tT+ «j  \ 


(64) 


J=> 


which  need  be  added  to  g  and  k  in  Eqs.  (54)  and  (55) 

e  e 


10 .  Bending  ( if  a  sj>l ieric.il  shell 

As  an  example  to  the  usefulness  of  the  she] 1  element  we  compute  the  large 
deformation  of  a  thin  spherical  shell  of  radius  1  and  v  =0.3  acted  upon 
by  a  surface  pressure  p  (negative  is  internal)  and  an  apex  force  !•'  (nega¬ 
tive  is  inward).  For  typographical  neatness  we  ignore  the  asterisks  on  p  and 
F.  The  geometry  and  loads  on  the  shell  are  shown  in  Fig.  3.  Because  of  sym¬ 
metry  we  consider  only  the  upper  half  of  the  sphere. 

Figure  1  shows  the  deformation  of  the  shell  discretized  with  14  elements, 
under  a  succession  of  polar  fortes  F  =  0,  -2.5,  -5.0,  ...,  -30  with  no  inter¬ 
nal  pressure.  Figure  2  shows  the  same  shell  under  the  same  loading  but  dis¬ 
cretized  with  24  elements. 


Tables  1  and  2,  referring  to  Figs.  1  and  2,  respectively,  list  Hgl^  as 
it  changes  with  the  Newton-Raphson  (NR)  iteration,  for  the  different  loadings. 
Evidently  the  computation  of  large  shell  deformations  can  become  expensive. 

Figure  3  shows  the  shell  of  Fig.  2  but  with  an  internal  pressure  p- 1 . C 
Figures  4  and  5  show  the  creation  of  a  dimple  in  a  thinner  (t  =  0.001)  shell 
The  solutions  for  F  =  10  and  F  =  20  in  Fig.  3  are  not  well  convergent. 
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.13E-9.46E0  .42E-2 

.  10E2 

.  35E1 

.46112 

.  44E2 

.  93E1 

. 68E 1  . 70E0 

. 17E4  . 16E0 

8 

. 45E- 1  .63E-5 

.  1  1E1 

.80112 

.  1 2  no 

.55111 

.  1  1 110 

.72114  .5211-1 

. 18E4  .86E-1 

9 

.3911-4  . 19E-9 

.  77E0 

.25110 

.52111 

.  14113 

.  1  3110 

.20113  .7611-4 

.22113  .26E-5 

10 

. 37E-9 

.  70 E- 2 

.  19112 

.6411-  3 

.  10111 

.2211-4 

.16113  .8211-9 

. 70E2  .25E-9 

11 

.  39  E- 4 

.13  E— 1 

.30E-5 

.  24E3 

. 48E-8 

.  33111 

.  22E1 

12 

.3411-9 

.6211-1 

.1211-9 

.  17111 

.3511-9 

.  33112 

.  5  1 E0 

13 

. 95E-6 

.  11111 

.  1  3110 

.1411-1 

14 

.2411-9 

.5411-1 

.27111 

.1711-4 

15 

.26F.-2 

.  46F.-4 

. 27E-9 

16 

.4211-6 

.  1  311-5 

17 

.2911-9 

.  1511-9 

I  Table  1:  Values  of  the  9  norm  of  the  global  gradient  HgH„  for  different 

t  1  z 

i 

j  i  appex  force  F,  as  they  diminish  with  the  Newt on-Raphson  (NR)  iteration.  This 


I  j  table  refers  to  Fig.  1. 


F 


-10.0  -12.5  -15.0  -17.5  -20.0  -22.5  -25.0  -27.50  -30. Q 


.  2 5 1-;  1  . 25E1  .  25K1  .  2 5 !•:  1  .  2 5 K 1  .251-1  .251-1  .25F.1  .25E1 

.151-4  .  32E4  .691-4  .  14E5  .21E5  .28E5  .461-5  .6  7E5  .82E5 

.  93E2  .  19 E 3  .45K3  .91E3  .97K3  .  1 3E4  .221-4  .21E4  .24E4 

.  53E2  .161-3  .  32  E  3  .43E3  .821-3  .  16E4  .221-4  .21E4  .44E4 

. 96E 1  .  14F.2  .  30 E 2  .45E2  .69E2  .  10E3  .  14E3  -  1 7E3  .  19E3 

.  12E2  . 26E2  .31E2  .29E2  .21E2  .49E2  . 11E3  .90E2  . 17E3 

.251-1  .  4  1 E 1  .94  El  .26E2  .451-2  .45E2  .  6  1  E2  .  10E3  .11E3 

.271.1  .  2 1E2  .  35 E2  .  13E2  .85E1  .15E2  .  1 2E2  .54E1  .63E1 

. 12EO  . 52EO  . 20E 1  .27E2  . 10E3  .84E2  . 18E3  . 78E3  .11E4 

'  .941-2  .36K 1  . 39K2  .46EI  .K8EO  .27E1  . 11  El  .24E1  .44E1 

.16K-5.10E-l.r3EO  .  38E2  .27K3  .45E3  .  14E4  .  11E4  .  15E4 

•44E-9.22E-2.36E1  .46E0  .91E0  .17E1  .66E1  .51E1  .61E1 

. 53E-8.12E-2  .23E2  .781-1  .  14E3  .211-2  .  12E3  .22E3 

.50F.-9  .34E-3  .I8E-I  .31K-1  .  18E0  .13E1  .22E0  .  33E0 

.5  IK-9  .  141-0  .101-0  .191-2  .  301-2  .181-3  .411-3 


13E2  .851-1  .  151-2  .  121-2  .54E1  .63E1 

27E2  .  10E3  .841-2  .  18E3  . 78E3  .11E4 

46 E 1  .881-0  .271-1  .111-1  .24E1  .44E1 

38E2  . 27E3  .45E3  . 14E4  . 11E4  . 15E4 
46E0  . 9 1E0  .  1 7E1  .66E1  . 5 1E1  .61E1 

23E2  .781-1  .  14E3  .211-2  .  12E3  .22E3 

181-1  .31E-1  . 18E0  . 1 3E1  .22E0  . 33EO 

141-0  .101-0  .191-2  .  301-2  .181-3  .411-3 

19E-5  .59E-5  .45E-2  .54E-1  .  20E0  .54E0 
13E-8  .33E-8  .31 E—  1  .45E1  . 35EO  .31E1 

38E-9  .34E -9  .65E-7  .41E-3  .65E-3  .18E-1 
.  29E-9  .56E-3  .11E-4.23E-1 
.44E-9  .73E-4  .11E-5 
. 40E-9 


Table  2.  Same  as  Table  1  but  Ne  =  24.  Refers  to  Fig.  2. 
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Figures 

1.  Thin  spherical  shell  bent  under  inwardly  directed  forces  F 
done  with  14  elements.  No  internal  pressure. 

2.  Same  as  Fig.  1  but  discretization  done  with  24  elements. 

3.  Same  as  Fig.  2  but  with  unit  internal  pressure. 

4.  A  dimple  in  a  thin  spherical  shell. 

3.  Same  as  Fig.  4  but  with  a  unit  internal  pressure. 


Discretization 


Ne=24 


Pmr»  in  Professional  Journal* 

Published  • 

1.  J.H.  Argyrls  and  I.  Fried:  The  LUMIMA  element  for  the  mstrix  displacement 
method.  J.  of  the  Royal  Aero.  Soc.,  514-517.  June  (1968). 

2.  J.H.  Argyrls,  1.  Fried  and  D.V.  Scharpf:  The  TET  20  and  TEA  8  element  for 
matrix  displacement  method.  J.  of  the  Royal  Aero.  Soc.  July  (1968). 

3.  J.H.  Argyrls,  I.  Fried  and  D.W.  Scharpf:  The  HERMES  8  element  for  the  matrix 
displacement  method.  J.  of  Royal  Aero.  Soc.,  613-617,  July  (1968). 

4.  1.  Fried:  Finite  element  analysis  of  time  dependent  phenomena.  AIAA  Journal 
2,  6,  (1969). 

5.  J.H.  Argyrls,  I.  Fried  and  D.V.  Scharpf:  The  TUBA  family  of  plate  elements 
for  the  matrix  displacement  method.  J.  of  the  Royal  Aero.  Soc.,  August  (1968). 

6.  I.  Fried:  More  on  gradient  iterative  methods  in  finite  element  analysis, 

AIAA  Journal,  7,  3,  March  (1969). 

7.  I.  Fried:  Gradient  iterative  methods  for  eigenproblems  in  the  finite  element 
analysis.  AIAA  Journal,  2*  4  (1969). 

8.  1.  Fried:  Some  aspects  of  the  natural  coordinate  system  in  the  finite  element 
method.  AIAA  Journal,  ]_■>  7,  July  (1969). 

9.  I.  Fried:  A  computational  procedure  for  the  solution  of  large  problems  arising 
from  the  finite  element  method.  Int.  J.  for  Numer.  Methods  in  Engineering,  2, 
*77-494  (1970). 

10.  X.  Fried:  Basic  computational  problems  in  the  finite  element  analysis  of 
shells.  Int.  J.  Solids  Structures,  Dec.  (1971). 

11.  I.  Fried:  Accuracy  of  finite  element  eigenproblems.  J.  Sound  and  Vibration, 

It  2.  (1971). 

12.  X.  Fried:  Discretisation  and  computational  errors  in  high  order  finite  elements. 
AIAA  Journal,  9,  10,  (1971). 

13.  I.  Fried:  N-step  conjugate  gradient  minimisation  scheme  for  non-quadratic 
functions.  AIAA  Journal,  2*  2286-2287  (1971). 


1*.  Z.  Fried:  Optimal  gradient  minimisation  scheme  for  finite  element  eigen- 
problems.  Journal  Sound  and  Vibration,  20,  3,  333-342  (1972). 


15.  X.  Fried:  Accuracy  of  complex  finite  elements.  AIAA  Journal.  10,  3, 

347-349  (1972). 

16.  Z.  Fried:  Condition  of  finite  element  matrices  generated  from  non-uniform 
mashes.  AIAA  Journal,  .10,  2,  219-221  (1972). 

17.  I.  Fried:  Bounds  on  the  extremal  eigenvalues  of  the  finite  element  matrices 
and  their  spectral  condition  number.  J.  Sound  and  Vibration,  22,  4,  407-418. 

18.  X.  Fried  and  Shok  Keng  Yang:  Best  finite  element  distribution  around  a 
singularity.  AIAA  Journal,  JJ),  3,  1244-1246  (1972). 

19.  X.  Fried:  Possible  loss  of  accuracy  In  curved  (isoparametric)  finite  elements. 
Journal  of  Sound  and  Vibration,  23,  4,  507-513  (1972). 

20.  1.  Fried:  Perturbation  errors  in  the  finite  element  method.  J.  Appl.  Mech. , 
629-631,  June  (1972). 

21.  X.  Fried:  Condensation  of  finite  element  eigenproblems.  AIAA  Journal,  10, 

11,  1529-1530,  November  (1972). 

22.  I.  Fried:  Shape  functions  and  the  accuracy  of  arch  finite  elements.  AIAA 
Journal,  U.,  3,  287-291  (1973). 

23.  X.  Fried:  Influence  of  Poisson's  ratio  on  the  condition  of  the  stiffness 
matrix.  Int.  J.  Solids  Structures,  9t  323-329  (1973). 

24.  1.  Fried:  Shear  in  C°  and  C*  bending  finite  elements.  Int.  J.  Solids  and 
Structures,  9,  449-460  (1973). 

25.  X.  Fried  and  Shok  Keng  Yang:  Triangular  nine-degrees-of-freemon,  C°,  plate 
bending  element  of  quadratic  accuracy.  Quart,  of  App.  Math.,  303-312,  Oct. 
(1973). 

26.  I*  Fried:  Finite  element  methods;  accuracy  at  a  point.  Quart,  of  App.  Math., 
149-161,  July  (1974). 

27.  I.  Fried:  Accuracy  and  condition  of  curved  (isoparametric)  finite  elements. 
Journal  of  Sound  and  Vibration,  31.*  345-357  (1973). 

28.  X.  Fried:  Boundary  and  interior  approximation  errors  in  the  finite  element 
method.  J.  App.  Mech.  Paper  No.  73-APM-A1A  (1973). 

29.  X.  Fried:  Bounds  on  the  spectral  and  maximum  norms  of  the  finite  element 
stiffness,  flexibility  and  mass  matrices.  Int.  J.  Solids  and  Structures,  9, 
1013-1034  (1973), 

30.  X.  Fried:  Residual  energy  balancing  technique  in  the  generation  of  plate  bending 
finite  elements.  Computers  and  Structures,  2,  771-778  (1974). 


31.  I.  Fried:  Numerical  integration  in  the  finite  element  method.  Computers 
and  Structures,  4,  921-932  (1974). 

32.  X.  Fried:  Note  on  the  finite  element  analysis  of  the  axisynmetric  elastic 
solid.  Int.  J.  Solids  and  Structures,  22*  383-386  (1974). 

33.  X.  Fried:  Finite  element  analysis  of  incompressible  material  by  residual 
energy  balancing.  Int.  J.  Solids  and  Structures,  J£,  993-1002  (1974). 

34.  I.  Fried  and  D.S.  Malkus:  Finite  element  mass  matrix  lumping  by  numerical 
integration  with  no  convergence  loss.  Int.  J.  Solids  and  Structures,  11. 
461-466  (1975). 

35.  X.  Fried:  Monoton  finite  element  matrices  and  their  computed  condition 
numbers.  Computers  and  Structures,  j>,  317-319  (1975). 

36.  I.  Fried:  Finite  element  analysis  of  thin  elastic  shells  with  residual 
energy  balancing  and  the  role  fo  rigid  body  modes.  J.  App.  Mech.  paper 
No.  75-AFM-6  (1975). 


37.  I.  Fried  and  J.  Metzler:  SOR  vs.  conjugate  gradients  in  a  finite  element 
discretization.  Xnt.  J.  Numer.  Methods  in  Eng.  12,  1329-1342  (1978). 

38.  X.  Fried  and  J.  Metzler:  Displacement,  strain  and  stress  error  nodal  lines 
in  finite  elements.  Computers  and  Structures,  9»,  335-339  (1978). 

39.  1.  Fried  and  J.  Metzler:  Conjugate  gradient  solution  of  a  finite  element 
elastic  problem  with  high  Poisson  ratio.  Comp.  Meth.  App.  Mech.,  Eng.  .15, 
1,  83-84  (1978). 


40. 

I.  Fried 

p-  »• 

41. 

it  Fried 
polutlon, 

42. 

I.  Fried 
626-628 

43. 

X.  Fried 
problems, 

44. 

X.  Fried 
illty  of 

45. 

I  Fried: 
estsnsib: 

Comp.  Meth.  App.  Mech. 


,  CMAME,  22,  229-240  (1980). 

nt  computations  of  the  equilibrium  and  stab- 

t  analysis  of 


46.  X.  Triad:  Stability  and  equilibrium  of  the  straight  and  curved  elastlea- 
flnlte  element  computation.  CMAMI,  28,  49-61  (1981). 


47.  I.  Fried:  Finite  element  computation  of  large  rubber  membrane  deformations. 
IJNME  18,  653-660  (1982). 

48.  I.  Fried:  Finite  element  computation  of  large  elastic  deformations:  Pro¬ 
ceedings  of  the  Brunei  Conference  of  the  mathematics  of  finite  elements 
and  applications,  MEFLAP  IV,  J.R.  Whiteman  Editor,  Academic  Press,  143-159 
(1982) . 

49.  I.  Fried:  Nonlinear  finite  element  computation  of  the  equilibrium  stability 
and  motion  of  the  extensional  beam  and  ring.  CMAME  38,  29-44  (1983). 

50.  I.  Fried:  Reflections  on  the  computational  approximation  of  elastic  in¬ 
compressibility.  Computers  and  Structures  17,  161-168  (1983). 


In  Print 

51.  I.  Fried:  Observations  on  stiffness  and  flexibility  modal  identification 
techniques.  CMAME  (1982). 

52.  I.  Fried:  Orthogonal  accession  to  the  equilibrium  curve.  CMAME  (1983). 

53.  I.  Fried:  On  conditionally  stable  explicit  time  integration  methods  in 
elastodynamics  and  heat  transfer.  CMAME  (1984). 

54.  I.  Fried:  Nonlinear  finite  element  analysis  of  the  thin  shell  of  revolu¬ 
tion.  CMAME  (1984). 


1.  Etude  de  flambage  d'une  coque  conique  sous  1'effet  combine  d'une  force  axiale 
•t  torsion  an  presence  d'une  presaion  interne;  Travail  elabore  dans  le  cadre 
d*un  program  d'obtention  de  dlplome  de  Maitre  es-science  Aeronautique 
(M.Sc.As.),  Ecole  National  Superleure  de  1' Aeronautique,  Paris,  Dec.  (1965). 

2.  Discretisation  and  round-off  errors  in  the  finite  element  analysis  of  el¬ 
liptic  boundary  value  problems;  Ph.D.  Massachusetts  Institute  of  Technology, 
June  (1971) . 

I.  Fried:  Numerical  Solution  of  Differential  Equations,  Academic  Press, 
tomtom,  New  York  (1979). 


Papers  In  Professional  Conferences 


1.  J.B.  Spooner,  O.E.  Bronlund,  I.  Fried,  J.H.  Argyris:  The  change  in  the  (two 
dimensional)  stress  distribution  round  a  rectangular  hole  with  rounded  cor¬ 
ners,  caused  by  a  varying  internal  pressure  and  temperature.  Proc.  XVII 
International  Aeronautical  Congress,  Belgrad,  24-29  Sept.  (1967). 

2.  I.  Fried:  A  computational  procedure  for  gradient  iterative  techniques  in 
the  finite  element  method.  EUROMECH  II  Conference,  Stuttgart,  8-10  April 
(1968). 

3.  J.H.  Argyris,  K.E.  Buck,  I.  Fried,  H.M.  Hilbert,  G.  Mareczek  and  D.W.  Scharpf: 
Some  new  elements  for  the  matrix  displacement  method.  2nd  Conference  on  Mat¬ 
rix  Methods  in  Structural  Mechanics,  Wright-Patterson  Air  Force  Base,  Ohio, 
15-17  October  (1968). 

4.  I.  Fried:  The  L2  and  L„  condition  numbers  of  the  finite  element  matrices 
and  the  pointwise  convergence  of  the  method.  Proc.  Conference  on  the 
Mathematics  of  the  Finite  Element  Method,  Edited  by  J.R.  Whiteman,  Aca¬ 
demic  Press,  Brunei  University,  18  April  (1972). 

5.  I.  Fried:  Finite  element  analysis  of  eigenproblems  -  theory  and  practice. 
National  Symposium  on  Computerized  Structural  Analysis  and  Design,  The 
George  Washington  University,  Washington,  D.C.  27-29,  March  (1972). 

6.  I.  Fried  and  D.S.  Malkus:  Energy  error  in  the  elastic  solution  when  an 
incompressible  solid  is  assumed  compressible.  Proceedings  of  the  U.S.- 
Germany  Seminar  on  the  Finite  Element  Method,  Edited  by  K..T.  Bathe, 

J.T.  Oden  and  W.  Wunderlich,  MIT  Press-,  August  (1976). 

7.  I.  Fried,  J.  Carson  and  Y.  Park:  Nonlinear  finite  element  analysis  of  gas 
flow  in  a  contact  bed  reactor.  Proc.  Conference  on  Finite  Element  Analysis 
in  Flow  Problems,  Portofino,  Italy,  16-18  June  (1976). 

8.  J.A.  Metzler  and  I.  Fried:  The  conjugate  gradient  method  with  finite 
elements.  2nd  IMACS  (AICA)  International  Symposium  on  Computer  Methods 
for  Partial  Differential  Equations,  Lehigh  University,  June  22-24  (1977). 

9.  I.  Fried  and  J.  Metzler:  Successive  overrelaxation  parameters  for  generrl 
finite  element  meshes,  1978  ACM  Computer  Science  Conference,  February  21-23, 
Detroit,  Michigan. 

10.  I.  Fried:  High  order  nonlinear  tangent  element  data  by  discrete  integration, 
Proc.  of  the  5th  Invitational  symposium  on  Finite  Elements,  Finite  Differ¬ 
ences  and  Calculus  of  Variations,  81-118,  H.  Kardestuncer  Editor,  University 
of  Connecticut,  May  2  (1980). 


*•■*,■»!*"* 


I 


AD- A 1 38  468  LINEAR  AND  NONLINEAR  FINITE  ELEHENTS(U)  BOSTON  UNIV  MA  V 
DEPT  OF  MATHEMATICS  I  FRIED  DEC  83  BU-1-84  J 

N000 14-76 -C- 0036 

UNCLASSIFIED  F/G  12/1  NL 


■SESSttL 

1*  X* Fried:  Finite  element  net bod  in  fluid  dynamics  end  heet  transfer.  Report 
■o.  38,  Znstltut  fur  stati.lt  und  Dynamik  der  Luft-und  Rauafhartkonstruktlonen, 
April  C1967 ) . 

2.  J.H.  Argyris,  W.  Boesard,  I.  Fried  and  H.M.  Hilbert:  A  fully  compatible 
plate  bending  element.  XSD  Report  No.  42,  Unlye rs It at  Stuttgart,  Deceaber 
(1967). 

3.  1*  Fried  and  D.S.  Malleus :  A  finite  element  displacement  model  valid  for  any 
value  of  the  compressibility.  Report  to  the  US  Office  of  Naval  Research 
(1976) . 

4.  A.  Johnson  and  I.  Fried:  High  order  nonlinear  finite  element  analysis  of 
the  axlsynmetric  rubber  membrane.  Report  6-80  to  the  Office  of  Naval  Re¬ 
search  (1980). 

Symposia 

1.  The  17th  International  Aeronautical  Congress.  Belgrad,  24-29  September, 

1967. 

2.  EUROMECH  11,  Experience  in  Computer-Based  Analysis  of  Three-Dimensional 
Media.  University  of  Stuttgart,  8-10  April  1968.  Invited. 

3.  Second  Conference  on  Matrix  Methods  in  Structural  Mechanics.  Wrlght- 
Patterson  Air  Force  Base,  Ohio,  15-17  October,  1968. 

4.  Seminar,  McGill  University,  Department  of  Electrical  Engineering.  Montreal, 
December  1971.  Invited. 

5.  Seminar  on  the  Role  of  Computers  In  Structural  Analysis,  Design  and  Op¬ 
timization.  Indian  Institute  of  Technology,  Madras,  December  11-13,  1972. 

6.  National  Symposium  on  Computerized  Structural  Analysis  and  Design.  George 
Washington  University,  March  27-29,  1972. 

7.  Boston  Numerical  Mathematics  Seminar.  Harvard  University,  April  12,  1972. 
Invited. 

8.  Seminar  on  the  Mathematical  Theory  of  Finite  Elements  Approximations. 

The  University  of  Alabama  in  Huntsville,  June  6-8,  1972.  Invited. 

9.  N aw  England  Bioengineering  Conference.  The  University  of  Vermont,  April  19-20, 
1973.  Invited. 

10.  US-Japan  Symposium  on  the  Mathematics  of  Finite  Elements  with  Emphasis  on 
Mom linear  Problems.  Lake  Hakone,  August,  1975.  Invited. 

11.  Applied  Mechanics  Conference.  Rensselaer  Polytechnic  Institute,  Troy, 

Seo  York,  June  23-25,  1975. 


.i’' 


12.  The  Mathematics  of  Finite  Elements  and  Applications.  Brunal  University , 
April  15-17.  1975. 

13.  National  Bureau  of  Standards  Mathematics  Colloquium  on  Deforsatlon  and 
Mechanical  Failures.  Gaithersburg.  May  16,  1975.  Invited. 

14.  IGAD-Intemational  Computer  Analysis  and  Design,  Second  Symposium  on  the 
Finite  Eleaient  Analysis  of  Fluid  Flow,  Santa  Margharlta  Ligure ,  Italy, 

June  14-18,  1975. 

15.  DS-Gexsany  Symposium  on  Finite  Element  Analysis.  Massachusetts  Institute 
of  Technology,  9-13  August  1976. 

16.  Second  IMACS  (AICA)  International  symposium  on  Computer  Methods  for  Partial 
Differential  Equations,  Lehigh  University,  June  22-24,  1977. 

17.  Fifth  Invitational  sympositm  on  Finite  Elements,  Finite  Differences  and  cal¬ 
culus  of  Variations,  University  of  Connecticut,  May  2  (1980).  Invited. 

18.  Fourth  conference  on  the  Mathematics  of  Finite  Elements  and  Applications, 
MEFLAP  IV,  Brunei  University,  May  1,  1981.  Invited. 

19.  Cambridge  symposium  on  optical  and  Electro-optical  Engineering,  Nov.  6-10 
(1983).  Invited. 
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Finite  elements,  Nonlinear  elasticity,  plate,  cable.  Rod, 
membrane,  shell  of  Revolution,  Incompressible,  Orthogonal 
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A  general  procedure  is*,  presented  for  the  setting  up  of  the  lle- 
ment  tangent  stiffness  matrix  and  element  gradient  through!  the 
discrete  sampling  of  the  nonquadratic  total  potential  enmrflBr* 
Computations' are  done  for  the  plate,  the  bent  rod,  the  cpUi 
the  rubber  membrane,  and  the  thin  shell  of  revolution. 
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